2.1. Wiener–Hopf equation

Inspired by Crighton (Reference Crighton1992), who investigated jet-edge tones, we focus on a two-dimensional jet flow of a vortex-sheet type in this paper. To render the model more realistic, the incident acoustic wave takes an asymptotic form generated through the interaction between shock and instability waves, as obtained previously by the authors (Li & Lyu Reference Li and Lyu2022, Reference Li and Lyu2023).

Schematic illustration is shown in figure 1. As we can see, a supersonic jet is issued from the nozzle and takes the form of a vortex sheet. The coordinate axes $(x,y)$ are chosen to be parallel and perpendicular to the nozzle centreline, respectively, while the origin point $O$ is located at the centre of the nozzle exit plane. According to Pack (Reference Pack1950), the base flow is chosen as the fully expanded jet velocity $\tilde {U}$ , which is different from the jet exit velocity $\tilde {U}_j$ . Here, $\tilde {D}$ is the height of the fully expanded jet. Note that $\tilde {D}$ is generally not equal to the height of the nozzle $\tilde {D}_0$ (Tam Reference Tam1972), but is less or greater than $\tilde {D}_0$ depending on whether the jet is overexpanded or underexpanded (Tam & Tanna Reference Tam and Tanna1982). The densities inside and outside the jet are denoted by $\tilde {\rho }_{\textrm {i}}$ and $\tilde {\rho }_{\textrm {o}}$ , respectively, while $\tilde {T}_{0,\textrm {i}}$ and $\tilde {T}_{\textrm {o}}$ represent the total reservoir temperature and the ambient temperature outside the jet, respectively. The nozzle lips, represented by two flat plates in figure 1, are assumed to be semi-infinite with a negligible thickness, extending to $-\infty$ in the $x$ direction. Downstream of the nozzle lip, an acoustic wave, with a speed of sound $\tilde {a}_\infty$ and a frequency of $\tilde {\omega }$ , is emitted from a source region and propagates upstream to the nozzle lip, triggering the instability waves.

To simplify the problem, we non-dimensionalise velocities and lengths by $\tilde {U}$ and $\tilde {D}$ , respectively. In addition, densities and temperatures are non-dimensionalised by $\tilde {\rho }_{\textrm {o}}$ and $\tilde {T}_{\textrm {o}}$ , respectively. Note that we use the symbols with a tilde to represent dimensional variables, while those without denote non-dimensional variables. For example, after the non-dimensionalisation, the jet exit velocity is $U_j$ and the base flow takes the form

(2.1) \begin{equation} \boldsymbol {u}_{0}= \begin{cases} 0, & |y|\gt 1/2, \\ {\boldsymbol{e}}_{x}, & |y|\leqslant 1/2, \end{cases} \end{equation}

where ${\boldsymbol{e}}_{x}$ is the unit vector in the $x$ direction.

Considering that the mean pressure outside and inside of the vortex sheet remains the same when a perfect gas is assumed, the density ratio between the outside and inside of the vortex sheet, here defined as $\sigma$ , is equal to $M_{\textrm {o}}^2 / M_{\textrm {i}}^2$ . Here, $M_{\textrm {o}}$ and $M_{\textrm {i}}$ are defined as $1/c_{\textrm {o}}$ and $1/c_{\textrm {i}}$ , respectively, where $c_{\textrm {o}}$ and $c_{\textrm {i}}$ represent the non-dimensional speed of sound outside and inside of the jet, respectively. The relation between the two Mach numbers is determined using Crocco–Busemann’s rule, i.e.

(2.2) \begin{equation} M_{\textrm {o}}=\dfrac {M_{\textrm {i}}}{\sqrt {1+\dfrac {\gamma -1}{2}M_{\textrm {i}}^2}}\nu ^{1/2}, \end{equation}

where $\nu$ denotes the temperature ratio between $\tilde {T}_{0,\textrm {i}}$ and $\tilde {T}_{\textrm {o}}$ , and $\gamma$ represents the specific heat ratio.

Figure 1.

Schematic of the vortex-sheet flow configuration and Cartesian coordinates. The origin is fixed at the centre of the nozzle while $x$ - and $y$ -axes represent the streamwise and cross-flow coordinates, respectively.

Before proceeding to solve the problem, it is useful to examine the designed transfer function using the Buckingham $\pi$ theorem. As mentioned above, eight physical variables are involved in this model, i.e. $\tilde {D}$ , $\tilde {U}$ , $\tilde {\rho }_{\textrm {i}}$ , $\tilde {\rho }_{\textrm {o}}$ , $\tilde {T}_{0,\textrm {i}}$ , $\tilde {T}_{\textrm {o}}$ , $\tilde {a}_\infty$ and the frequency of the acoustic wave $\tilde {\omega }$ , while four independent dimensions can be found. Employing the Buckingham $\pi$ theorem, it is evident that the transfer function $H$ between the newly excited instability wave and the acoustic waves is related to four non-dimensional parameters. These four parameters are chosen to be the Mach number $M_{\textrm {i}}$ , the Strouhal number $\omega$ , the temperature ratio $\nu$ and the density ratio $\sigma$ . The density ratio can be calculated using $\nu$ and $M_{\textrm {i}}$ through the Crocco–Busemann rule (2.2). Therefore, there remain three independent non-dimensional parameters and our objective is to find the transfer function $H (M_{\textrm {i}},\omega ,\nu )$ in the subsequent analysis.

To have a realistic incident acoustic wave, we take the result obtained in the SII work (Li & Lyu Reference Li and Lyu2022, Reference Li and Lyu2023). Pack’s model was used to describe shock waves, while instability waves were calculated using the linear stability analysis. The resulting non-dimensionalised velocity potential function of the acoustic waves due to one-cell interaction can be obtained by performing the inverse Fourier transform, i.e.

(2.3) \begin{equation} \phi _{\textit{in}}(x,y)=\frac {1}{2\pi }\int _{-\infty }^{\infty }D_{\textrm {o}}(\lambda )\textrm {exp}\left(-\left(\textrm {i}\lambda x+\sqrt {\lambda ^2-\omega ^2 M_{\textrm {o}}^2}y\right)\right)\textrm {d}\lambda . \end{equation}

Herein, $\lambda$ denotes the streamwise wavenumber and $D_{\textrm {o}}(\lambda )$ represents a coefficient, which can be readily calculated when $M_{\textrm {o}}$ and $\omega$ are provided. The detailed expression of $D_{\textrm {o}}(\lambda )$ can be found in Appendix C. The integration in (2.3) can be estimated in the far field using the saddle point method. In this far field, the primary contribution to the integral value comes from the vicinity of the saddle point $\lambda _0=-\omega M_{\textrm {o}}\textrm {cos} \theta$ . As a result, the potential function can be written as follows:

(2.4) \begin{equation} \phi _{\textit{in}}(r,\theta )=\sqrt {\frac {{\omega M_{\textrm {o}}}}{{2\pi }}}\frac {\textrm {e}^{\textrm {i}\omega (M_{\textrm {o}}r-t-\pi /4)}}{\sqrt {r}}D_{\textrm {o}}(-\omega M_{\textrm {o}}\textrm {cos} \theta )\textrm {sin} \theta , \end{equation}

where $\theta$ denotes the observer angle with respect to the downstream direction and $r$ represents the distance from the screech source region to the nozzle lip. Further details on the saddle point method can be found from Crighton et al. (Reference Crighton, Dowling, Williams, Heckl and Leppington1992).

As shown in (2.4), $\phi _{ {in}}(r,\theta )$ depends on $r$ . Experimental observations indicated that the screech source was situated several shock cells downstream from the nozzle exit (Suda et al. Reference Suda, Manning and Kaji1993; Kaji & Nishijima Reference Kaji and Nishijima1996; Malla & Gutmark Reference Malla and Gutmark2017). Therefore, the streamwise distance between the source position to the nozzle lip is written as

(2.5) \begin{equation} h=n s, \end{equation}

where $n$ typically takes the value of 3, 4 or 5 (Suda et al. Reference Suda, Manning and Kaji1993; Li & Gao Reference Li and Gao2010; Malla & Gutmark Reference Malla and Gutmark2017) and $s$ denotes the shock spacing, which can be calculated using Pack’s model (Pack Reference Pack1950; Tam Reference Tam1986). Note that we assume the shock cells are equally spaced downstream of the jet flow, while in realistic jets, the shock spacing gradually decreases (Norum & Seiner Reference Norum and Seiner1982). This non-uniform spacing may change the screech frequency prediction slightly. However, in this paper, our focus is to develop a transfer function at a given frequency. Such approximation should suffice for such a purpose.

We consider the potential function of the upstream forcing $\Phi _{ {in}}(r,\theta )$ in the proximity of the nozzle lip, where the value of $y-1/2$ is relatively small compared with $|x|=h$ . Therefore, the distance $r$ in (2.4) can be approximated as $r\approx h$ . Similarly, we have $\sin \theta \approx (y-1/2)/h$ and $\cos \theta \approx 1$ . Consequently, in the vicinity of the nozzle, the potential function $\Phi _{ {in}}$ can be reformulated to a similar form to that of Crighton (Reference Crighton1992), i.e.

(2.6) \begin{equation} \Phi _{\textit{in}}(Y)=GY, \end{equation}

where

(2.7) \begin{equation} G=\sqrt {\frac {{\omega M_{\textrm {o}}}}{{2\pi }h^{3}}}D_{\textrm {o}}(\omega M_{\textrm {o}}) \textrm {e}^{\textrm {i}\omega (M_{\textrm {o}}h-\pi /4)}, \end{equation}

(2.8) \begin{equation} Y=y-\tfrac{1}{2}. \end{equation}

The goal is therefore to find the response of the jet flow due to the upstream forcing $GY$ near the nozzle lip. In this paper, we focus on the upper half-plane ( $y\gt 0$ ), while a similar analysis can be carried out for the lower half-plane. Considering that the time-harmonic assumption is used, the factor $\textrm {e}^{-\textrm {i}\omega t}$ in all the time-dependent fields is omitted for clarity in the rest of this paper. Because we assume the initial base flow to be irrotational and inviscid both inside and outside of the vortex sheet (Batchelor & Gill Reference Batchelor and Gill1962), the velocity potential induced by the incident acoustic wave near the nozzle lip can be described by a potential function. Therefore, taking into account (2.6), the total velocity potentials outside and inside the jet flow can be defined by

(2.9) \begin{equation} \Phi =GY+\phi (x,y), \end{equation}

(2.10) \begin{equation} \Psi =GY+\psi (x,y), \end{equation}

respectively, where $\phi$ and $\psi$ are the corresponding potential functions of the scattered sound fields in the region $y\gt 1/2$ and $0\leqslant y\leq 1/2$ , respectively.

It is known that $\phi$ and $\psi$ both satisfy the convective wave equation, i.e.

(2.11) \begin{equation} {\nabla }^{2}\phi +\omega ^2 M_{\textrm {o}}^2\phi =0, \end{equation}

(2.12) \begin{equation} {\nabla }^{2}\psi +\left ( \omega + \textrm {i}M_{\textrm {i}}\frac {\partial }{\partial x}\right )^2\psi =0. \end{equation}

Hard-wall conditions are imposed on $y=1/2$ for $-\infty \lt x \leq 0$ . Note that, strictly speaking, these conditions should be applied at $y=\tilde {D_0}/2\tilde {D}$ , which can be either less than or greater than $1/2$ depending on whether the jet is underexpanded or overexpanded (Tam & Tanna Reference Tam and Tanna1982). However, as mentioned earlier, the difference between $\tilde {D_0}/2\tilde {D}$ and $1/2$ is relatively small compared with the distance $r$ . To enable analysis progress, we neglect this difference and assume that the plate is nearly parallel to the boundary of the jet. Therefore, we obtain

(2.13) \begin{equation} \frac {\partial \phi }{\partial y}=\frac {\partial \psi }{\partial y}=-G \quad (y=1/2,\ -\infty \lt x\leq 0). \end{equation}

The absence of decay in the forcing term $-G$ leads to convergence issues when the Fourier transform is applied to (2.13). To circumvent this problem, Crighton (Reference Crighton1992) replaced (2.13) by

(2.14) \begin{equation} \frac {\partial \phi }{\partial y}=\frac {\partial \psi }{\partial y}=-G\textrm {e}^{\epsilon x}, \end{equation}

where $\epsilon$ is a real and positive term, which retains the physical characteristics of the upstream forcing near the nozzle lip but removes the forcing further upstream. We follow the same treatment here. Note that in realistic jets, the acoustic wave decays algebraically as it propagates from the downstream source region to the upstream nozzle lip. In this study, we model the behaviour of the acoustic wave near the nozzle lip using the exponential decay rate, which will be uniquely determined in § 2.3.

The kinematic and dynamic boundary conditions can be used on the vortex sheet at $y=1/2,\ 0\lt x\lt \infty$ , i.e.

(2.15) \begin{equation} \left (-\textrm {i} \omega +\frac {\partial }{\partial x}\right )\frac {\partial \phi }{\partial y}=-\textrm {i}\omega \frac {\partial \psi }{\partial y}, \end{equation}

(2.16) \begin{equation} -\textrm {i}\omega \sigma \phi = \left (-\textrm {i} \omega +\frac {\partial }{\partial x}\right )\psi . \end{equation}

Fourier transforms $\phi (k,y), \psi (k,y)$ are defined by

(2.17) \begin{gather} \phi (k,y)=\int _{-\infty }^{+\infty } \phi (x,y)\textrm {e}^{\textrm {i} k x}\textrm {d} x, \end{gather}

(2.18) \begin{gather} \psi (k,y)=\int _{-\infty }^{+\infty } \psi (x,y)\textrm {e}^{\textrm {i} k x}\textrm {d} x, \end{gather}

where $k$ is the wavenumber in the streamwise direction. Multiplying (2.11) and (2.12) with $\textrm {e}^{\textrm {i}kx}$ and subsequently integrating both expressions with respect to the variable $x$ , solutions to (2.11) and (2.12) can be found upon invoking the far-field boundary condition, i.e.

(2.19) \begin{equation} \phi (k,y)=\phi ^{*}(k)\textrm {e}^{-\gamma _1(k)y}, \end{equation}

(2.20) \begin{equation} \psi (k,y)=\psi _{\textrm {s}}^{*}(k)\sinh \gamma _2(k)y+\psi _{\textrm {c}}^{*}(k)\cosh \gamma _2(k)y, \end{equation}

where

(2.21) \begin{equation} \gamma _2(k)=\sqrt {k^2-(\omega +k)^2M_{\textrm {i}}^2} \quad \textrm {and} \quad \gamma _1(k)=\sqrt {k^2-\omega ^2 M_{\textrm {o}}^2}. \end{equation}

Herein, $\phi ^*(k)$ , $\psi _{\textrm {s}}^*(k)$ and $\psi _{\textrm {c}}^{*}(k)$ are functions of $k$ , whose detailed forms will be specified subsequently. To satisfy the far-field condition, the real part of $\gamma _{1}$ should be positive and its branch cut will be specified later. It is well known that both symmetric (varicose) and antisymmetric (sinuous) instability modes can exist inside the jet. However, for jet flows from high-aspect-ratio rectangular nozzles, the sinuous mode is dominant (Edgington-Mitchell Reference Edgington-Mitchell2019). Therefore, we consider the antisymmetric mode in this study, while the symmetric mode can be examined in a similar way should interest arise. In this case, $\psi (k,y)$ reduces to

(2.22) \begin{equation} \psi (k,y)=\psi ^{*}_{\textrm {s}}(k)\sinh \gamma _2(k)y. \end{equation}

In what follows, we use the Wiener–Hopf method to solve this boundary-value problem specified by (2.19), (2.22), and boundary conditions (2.14), (2.15) and (2.16). The details of the Wiener–Hopf technique can be found from Noble (Reference Noble1958). Using (2.14), we can obtain

(2.23) \begin{equation} \phi ^{\prime }(k,1/2)=\phi _{+}^{\prime }(k,1/2)+\frac {\textrm {i}G}{k-\textrm {i}\epsilon }, \end{equation}

where the prime represents the first derivative with respect to $y$ and $\phi _{+}(k,1/2)$ is defined by

(2.24) \begin{equation} \phi _{+}(k,1/2)=\int _{0}^{\infty }\phi (x,1/2)\textrm {e}^{\textrm {i}k x}\textrm {d}x. \end{equation}

Considering the kinematic boundary condition across the vortex sheet, i.e. (2.15), and the boundary conditions imposed on the nozzle lip, i.e. (2.14), we have

(2.25) \begin{equation} (k-\textrm {i}\epsilon )\psi ^{*}_{\textrm {s}}(k)\cosh \gamma _2(k)/2-\frac {G}{\omega }\epsilon =-\left (1+\frac {k}{\omega }\right )(k-\textrm {i}\epsilon )\frac {\gamma _1(k)}{\gamma _2(k)}\phi ^{*}(k)\textrm {e}^{-\gamma _1(k)/2}. \end{equation}

Following Crighton (Reference Crighton1992), the term $\dfrac {G}{\omega }\epsilon$ on the left-hand side of (2.25) is omitted, considering that $\epsilon$ is a small term.

Now, we define a new function related to the difference between the pressures outside and inside the vortex sheet as $D(k,1/2)$ , i.e.

(2.26) \begin{equation} D(k,1/2)=\int _{-\infty }^{\infty } \left (-\textrm {i}\omega \sigma \phi (x,1/2)-\left (-\textrm {i}\omega +\frac {\partial }{\partial x}\right )\psi (x,1/2)\right )\textrm {e}^{\textrm {i}k x}\textrm {d}x. \end{equation}

Due to the continuity of the pressure across the vortex sheet, i.e. (2.16), the pressure difference vanishes when $x\gt 0$ . Therefore, $D(k,1/2)$ reduces to

(2.27) \begin{equation} \begin{aligned} D(k,1/2) &= -\textrm {i}\omega \sigma \phi _{-}(k,1/2)-(-\textrm {i}\omega -\textrm {i}k)\psi _{-}(k,1/2)+\psi _0, \end{aligned} \end{equation}

where

(2.28) \begin{gather} \phi _{-}(k,1/2)=\int _{-\infty }^{0}\phi (x,1/2)\textrm {e}^{\textrm {i}k x}\textrm {d}x \quad \textrm {and} \quad \psi _{-}(k,1/2)=\int _{-\infty }^{0}\psi (x,1/2)\textrm {e}^{\textrm {i}k x}\textrm {d}x, \end{gather}

and $\psi _0$ denotes the finite value of $\psi (x,y)$ at the position $(0,1/2)$ . Substituting (2.19) and (2.22) into (2.27), we obtain

(2.29) \begin{equation} D(k,1/2)=-\textrm {i}\omega \sigma \phi ^{*}(k)\textrm {e}^{-\gamma _1(k)/2}-(-\textrm {i}\omega -\textrm {i}k)\psi ^{*}_{\textrm {s}}(k)\sinh \gamma _2(k)/2. \end{equation}

Combining (2.23), (2.25), (2.27) and (2.29) yields

(2.30) \begin{equation} \phi _{+}^{\prime }(k,1/2)+\frac {\textrm {i}G}{k-\textrm {i}\epsilon }=K(k)F_{-}(k), \end{equation}

where

(2.31) \begin{equation} K(k)=\frac {1}{\dfrac {\sigma }{\gamma _1(k)}+\dfrac {(1+k/\omega )^2\tanh \gamma _2(k)/2}{\gamma _2(k)}}, \end{equation}

(2.32) \begin{equation} F_{-}(k)=(1+k/{\omega })\psi _{-}(k,1/2)-\sigma \phi _{-}(k,1/2)-{\textrm {i}}\psi _0/{\omega }. \end{equation}

Equation (2.30) represents the Wiener–Hopf equation we aim to solve, with the kernel given by $K(k)$ . Here, $1/K(k)=0$ represents the dispersion relation of the antisymmetric mode of the jet instability wave. It can be seen that as $|k|\rightarrow \infty$ , $|K(k)|\rightarrow |k|^{-1}$ . Now, following the procedures of the Wiener–Hopf method, we assume that the kernel can be decomposed into

(2.33) \begin{equation} K(k)=K_{+}(k)K_{-}(k), \end{equation}

where $K_{+}(k)$ and $K_{-}(k)$ are analytic and non-zero in the upper and lower half-planes of $k$ , respectively. Meanwhile, both $K_{+}(k)$ and $K_{-}(k)$ behave as $O(k^{-1/2})$ as $k$ approaches infinity. The actual decomposition will be shown in § 2.2. Routine use of Wiener–Hopf decomposition leads to

(2.34) \begin{equation} \frac {\phi _{+}^{\prime }(k,1/2)}{K_{+}(k)}+\left (\frac {1}{K_{+}(k)}-\frac {1}{K_{+}(\textrm {i}\epsilon )}\right )\frac {\textrm {i}G}{k-\textrm {i}\epsilon }=E(k), \end{equation}

(2.35) \begin{equation} K_{-}(k)F_{-}(k)-\frac {1}{K_{+}(\textrm {i}\epsilon )}\frac {\textrm {i}G}{k-\textrm {i}\epsilon }=E(k). \end{equation}

Here, $E(k)$ is an entire function of $k$ . Due to the Kutta condition (Crighton Reference Crighton1985), $\phi (x,1/2)$ behaves like $x^{3/2}$ near the nozzle lip. The corresponding value of $ \phi ^\prime (x,1/2)$ is therefore $O(x^{1/2})$ . It follows that $\phi _{+}^{\prime }(k,1/2)\rightarrow k^{-3/2}$ as $k\rightarrow \infty$ (Erfelyi Reference Erfelyi1958). Thus, from Liouville’s theorem,

(2.36) \begin{equation} E(k)\equiv 0. \end{equation}

From (2.23) and (2.34), $\phi ^{\prime }(k,1/2)$ can be readily solved, i.e.

(2.37) \begin{equation} \phi ^{\prime }(k,1/2)=\frac {\textrm {i}G}{k-\textrm {i}\epsilon }\frac {K_{+}(k)}{K_{+}(\textrm {i}\epsilon )}. \end{equation}

It can be seen that $\phi ^{\prime }(k,1/2)$ is $O(k^{-3/2})$ as $|k|$ approaches infinity in the upper half-plane. The corresponding velocity component $\phi ^\prime (x,y)$ behaves as $x^{1/2}$ when $x\rightarrow 0$ , which leads to singularity in $\psi ^\prime (x,y)$ due to (2.15). Therefore, (2.37) leads to singular behaviour when $|k|\rightarrow \infty$ , where (2.37) can be approximated by (Crighton Reference Crighton1992)

(2.38) \begin{equation} \frac {\textrm {i}G}{K_{+}(\textrm {i}\epsilon )}\frac {K_{+}(k)}{k}. \end{equation}

This part will be further discussed in § 2.3, where the unsteady Kutta condition is used to determine the new instability wave by removing the emerging singularity.

2.2. Kernel decomposition

In this section, we factor the kernel $K(k)$ into $K_{+}(k)$ and $K_{-}(k)$ , which is the essential step of the Wiener–Hopf technique. Before proceeding with the decomposition of the kernel, it is necessary to properly define the branch cuts of $K(k)$ . To ensure that the real part of $\gamma _{1}$ is positive as $|k|\rightarrow \infty$ along the integration path (which becomes relevant when conducting the inverse Fourier transform if we are interested in calculating the scattered sound field), the branch cuts of $\gamma _{1}$ passing the branch points $k=\pm \omega M_{\textrm {o}}$ are chosen to extend to the upper and lower half-plane, respectively, as shown in figure 2(a). The branch is selected in such a way that $\textrm {arg}(\gamma _1)=-\pi /2$ at $k=\textrm {i}$ . However, the branch points of $\gamma _{2}$ are determined as $M_{1}=-\omega M_{\textrm {i}}/(M_{\textrm {i}}+1)$ and $M_{2}=-\omega M_{\textrm {i}}/(M_{\textrm {i}}-1)$ . The branch cuts passing $M_{1}$ and $M_{2}$ are chosen to extend to the lower and upper half-plane, respectively, for the convenience of the kernel decomposition. The branch of $\gamma _2$ is selected to ensure that $\textrm {arg}(\gamma _2)=-\pi /2$ at $k=\textrm {i}$ . Due to the choice of the branch cuts, we can find that $\gamma _2(k)/\gamma _1(k)\rightarrow \sqrt {1-M_{\textrm {i}}^2}$ , as $\omega \rightarrow 0$ and $k\rightarrow O(1)$ .

Figure 2.

(a) Schematic of the branch points, branch cuts and the integration path in the complex $k$ plane. (b) Schematic of the deformed integration path in the complex $t$ plane.

The kernel can be rewritten as

(2.39) \begin{equation} K(k)=\frac {\gamma _2(k)}{\sigma \gamma _2(k)/\gamma _1(k)+(1+k/\omega )^2\tanh \gamma _2(k)/2}. \end{equation}

A proper decomposition of $\gamma _2(k)$ is evident, i.e. $\gamma _2(k)=\sqrt {k-M_1}\sqrt {k-M_2}$ . Thus, in what follows, we only focus on the denominator, which we refer to as $Q(k)$ . It can be seen that $Q(k)$ depends on the dispersion relation of the antisymmetric mode of the jet instability wave. This dispersion relation describes several types of waves present in the jet, determined by its structural characteristics such as branches, zeros and poles. These information can be maintained if the kernel is decomposed by means of matched asymptotic expansions (Crighton Reference Crighton1992, Reference Crighton2001). To use this strategy, a small parameter should be specified as the basis for the expansion. This parameter should not only be small but can also be contained by the kernel. We introduce $S=\omega (M_{\textrm {i}}-1)^2$ as this parameter. In what follows, the complicated kernel is replaced by a simpler one by asymptotic expansions based on $S$ . Overlapping approximations are subsequently applied to guarantee the validity of the decomposition over the whole range of the wavenumber $k$ . The final result is obtained by a multiplicative composite approximation.

We first assume $k=O(1)$ and let $\omega \rightarrow 0$ . The first term of $Q(k)$ , i.e. $\sigma \gamma _2(k)/\gamma _1(k)$ , reduces to $\sigma \sqrt {1-M_{\textrm {i}}^2}$ and $Q(k)$ can be approximated as

(2.40) \begin{equation} \begin{aligned} Q(k) &\approx \left (\frac {k}{\omega }\right )^2\tanh \sqrt {\left(1-M_{\textrm {i}}^2\right)k^2}/2. \end{aligned} \end{equation}

The term $\sqrt {(1-M_{\textrm {i}}^2)k^2}$ is supposed to maintain the branch cut chosen for $\gamma _2(k)$ . For brevity, we rewrite $\sqrt {(1-M_{\textrm {i}}^2)k^2}$ as $\sqrt {1-M_{\textrm {i}}^2k^{*}}$ , where the script $^*$ represents that this parameter retains the specified branch cuts.

Now, we introduce the small parameter $S=\omega (M_{\textrm {i}}-1)^2$ and a new scaled wavenumber according to

(2.41) \begin{equation} q=k/S^{\frac {3}{4}}. \end{equation}

By assuming $S\rightarrow 0$ and $q\rightarrow O(1)$ , we can expand $Q(k)$ in (2.39) based on the small parameter $S$ (using methods of Taylor expansion), i.e.

(2.42) \begin{equation} \begin{aligned} Q(k)&\approx \sqrt {1-M_{\textrm {i}}^2}\left(\sigma -\frac {\sigma M_{\textrm {i}}^2S^{1/4}}{(M_{\textrm {i}}-1)^2\left(1-M_{\textrm {i}}^2\right)q^{*}}+\frac {1}{2}(M_{\textrm {i}}-1)^4S^{1/4}q^{2}q^{*}\right), \end{aligned} \end{equation}

where $q^*=k^*/S^{\frac {3}{4}}$ . Note that despite the small parameter $S\sim o(1)$ , the quantities $S^{1/4}\sim O(1)$ . To put this into perspective, let us consider an example with $M_{\textrm {i}}=1.1$ . At the screech frequency, employing Powell’s original model yields $S=0.026$ (with $\kappa =0.6$ and $\nu =1$ ). However, $S^{1/4}=0.40$ , which is of the order of unity. The scale separation under other jet Mach numbers is discussed in Appendix A. This demonstrates that although $S$ itself may be small, its corresponding fractional power can still have magnitudes of the order of unity in our interested application.

We can show that (2.40) and (2.42) do overlap, as $k\rightarrow 0$ and $q\rightarrow \infty$ , respectively, with common value $({1}/{2})\sqrt {1-M_{\textrm {i}}^2}(M_{\textrm {i}}-1)^4S^{1/4}q^{2}q^{*}$ . A multiplicative composite approximation then reads

(2.43) \begin{equation} \begin{aligned} Q(k)&\approx \sqrt {1-M_{\textrm {i}}^2} \frac {\tanh \dfrac {1}{2}\sqrt {1-M_{\textrm {i}}^2} k}{\dfrac {1}{2}\sqrt {1-M_{\textrm {i}}^2} k} \frac {1}{q^{*}}R(q), \end{aligned} \end{equation}

where

(2.44) \begin{equation} R(q)=\sigma q^{*}-\frac {\sigma M_{\textrm {i}}^2S^{1/4}}{(M_{\textrm {i}}-1)^2\left(1-M_{\textrm {i}}^2\right)}+\frac {1}{2}(M_{\textrm {i}}-1)^4S^{1/4}q^{4}. \end{equation}

This approximation strategy enables direct calculation of the wavenumber for the sinuous mode of the jet instability wave. To validate the decomposition, the wavenumber analytically obtained by letting $Q(k)$ in (2.43) equal to 0 is compared with that numerically solved from the dispersion relation $1/K(k) = 0$ . As shown in figure 3, the agreement between the two wavenumbers is good when $M_{\textrm {i}} \gt 1.2$ , which corresponds to the screech frequency in the range of $\omega \lesssim 1.8$ . This indicates that to accurately estimate the wavenumber, frequency should be a small term. Additional validations under other operation conditions are provided in Appendix A.

Figure 3.

Comparison between the wavenumbers of the sinuous mode of the jet instability wave obtained by the decomposed kernel (2.43) and the dispersion relation. $\alpha ^{*}$ and $\alpha ^{v}$ denote the wavenumber calculated from the decomposed kernel and $M(k)$ , respectively, while the subscripts $_r$ and $_i$ represent the real and imaginary part of the parameter, respectively. $\sigma$ is assumed to be $1$ .

We are now in a position to factor the function $Q(k)$ . From Noble (Reference Noble1958), $({\tanh \sqrt {1-M_{\textrm {i}}^2} k/2})/({\sqrt {1-M_{\textrm {i}}^2} k/2})$ can be decomposed to

(2.45) \begin{equation} \begin{aligned} \frac {\tanh \dfrac {1}{2}\sqrt {1-M_{\textrm {i}}^2} k}{\dfrac {1}{2}\sqrt {1-M_{\textrm {i}}^2} k} &= T_{+}(k)T_{-}(k), \end{aligned} \end{equation}

where

(2.46) \begin{gather} T_+(k)=\frac {\Gamma \left(1/2+\sqrt {M_{\textrm {i}}^2-1}k/2\pi \right)}{\pi ^{1/2}\Gamma \left(1+\sqrt {M_{\textrm {i}}^2-1}k/2\pi \right)},\quad T_-(k)=\frac {\Gamma \left(1/2-\sqrt {M_{\textrm {i}}^2-1}k/2\pi \right)}{\pi ^{1/2}\Gamma \left(1-\sqrt {M_{\textrm {i}}^2-1}k/2\pi \right)}. \end{gather}

Here, $\Gamma$ represents the Gamma function. The decomposition of $1/q^*$ reads

(2.47) \begin{equation} \begin{aligned} \frac {1}{q^*} &=\mathcal {Q}_{+}(q)\mathcal {Q}_{-}(q), \end{aligned} \end{equation}

where

(2.48) \begin{equation} \mathcal {Q}_{+}(q)=\frac {1}{\sqrt {q+0}},\quad \mathcal {Q}_{-}(q)=\frac {1}{\sqrt {q-0}}. \end{equation}

Naturally, the branch cut of $\mathcal {Q}_{+}$ and $\mathcal {Q}_{-}$ is from $-0$ to $-\textrm {i}\infty$ and from $+0$ to $+\textrm {i}\infty$ , respectively.

The decomposition of $R(q)$ , however, is very difficult due to the presence of $q^{*}$ in (2.44).

To decompose $R(q)$ , we let $q$ be real first (Crighton Reference Crighton1992), so that $q^{*}=|q|$ . This is reasonable considering that the overlap of the two half-planes will converge to the real axis in the complex $k$ plane. Assuming that $R(q)$ can be decomposed to $R(q)=R_{+}(q)R_{-}(q)$ and taking the first derivative of $\textrm {ln}R(q)$ , we obtain

(2.49) \begin{equation} \frac {\textrm {d}}{\textrm {d}q}\textrm {ln}R_{+}(q)+\frac {\textrm {d}}{\textrm {d}q}\textrm {ln}R_{-}(q)=\frac {4q^3+a_1 |q|/q}{q^4+a_1|q|+a_0}, \end{equation}

where

(2.50) \begin{align} a_0=2\sigma M_{\textrm {i}}^2(M_{\textrm {i}}-1)^{-6}\left(M_{\textrm {i}}^2-1\right)^{-1},\nonumber \\ a_1=2\sigma (M_{\textrm {i}}-1)^{-4}S^{-1/4}.\qquad \end{align}

Multiplying the numerator and the denominator of the right-hand side term in (2.49) by $q^4-a_1|q|+a_0$ , (2.49) can be written as

(2.51) \begin{equation} \frac {\textrm {d}}{\textrm {d}q}\textrm {ln}R_{+}(q)+\frac {\textrm {d}}{\textrm {d}q}\textrm {ln}R_{-}(q)=\frac {(4q^3+a_1 |q|/q)(q^4-a_1|q|+a_0)}{(q^4+a_0)^2-a_1^2q^2}. \end{equation}

The zeros of the denominator $(q^4+a_0)^2-a_1^2q^2$ are defined as $q_n, n=1,\ldots ,8$ , which can be obtained by Ferrari’s formula. The details are provided in Appendix A. With the knowledge of these zeros, we can now reorganise (2.51) as

(2.52) \begin{equation} \begin{aligned} \frac {\textrm {d}}{\textrm {d} q} \ln \left (\frac {R_{+}(q)}{\prod _{n=5}^8\left (q-q_n\right )^{1/2}}\right )+\frac {\textrm {d}}{\textrm {d} q} \ln & \left (\frac {R_{-}(q)}{\prod _{n=1}^4\left (q-q_n\right )^{1/2}}\right ) \\ &=\sum _{n=1}^4 \xi _n \frac {|q|}{\left (q-q_n\right )}+ \sum _{n=5}^8\xi _n \frac {|q|}{\left (q-q_n\right )}+\frac {a_1 }{a_0 }\frac {|q|}{q}, \end{aligned} \end{equation}

where

(2.53) \begin{equation} \xi _n =\frac {-3a_1q_n^4+a_0a_1}{2q_n^2\left(4q_n^6+4a_0q_n^2-a_1^2\right)}. \end{equation}

The detailed derivation from (2.51) to (2.52) can be found in Appendix B. For brevity, we can define $\xi _9=a_1/a_0$ and $q_9=0$ . With these definitions, it is straightforward to find

(2.54) \begin{equation} \sum _{n=1}^9 \xi _n q_n=\sum _{n=1}^9 \xi _n=0. \end{equation}

An addictive decomposition of $|q|$ is now introduced, i.e.

(2.55) \begin{eqnarray} |q|&=&Q_{+}(q)+Q_{-}(q), \end{eqnarray}

where

(2.56) \begin{equation} Q_{+}(q)=q/2+\frac {\textrm {i}q}{\pi }\ln _{+}q,\quad Q_{-}(q)=q/2-\frac {\textrm {i}q}{\pi }\ln _{-}q. \end{equation}

The branch cuts of $\ln _{+}$ and $\ln _{-}$ extend from $0$ to $-\textrm {i}\infty$ and from $0$ to $+\textrm {i}\infty$ , respectively. The choice of the branch cut leads to

(2.57) \begin{equation} Q_+(q)=Q_{-}(-q). \end{equation}

Using the addictive decomposition (2.55) and the two identities (2.54) and (2.57), the right-hand side of (2.52) can obtain a similar addictive factorisation and $R_{+}(q)$ can be finally found as

(2.58) \begin{equation} R_{+}(q)=\prod _{n=5}^8(q-q_n)^{1/2}\textrm {exp}\left (-\frac {\textrm i}{\pi }\sum _{n=1}^{8}\xi _n q_n\int _{q/q_n }^{\infty }\frac {\ln _{+}t}{t-1}\textrm {d}t +\sum _{n=5}^{8}{\xi _n q_n}\int _{q/q_n }^{\infty }\frac {1}{t-1}\textrm {d}t \right ). \end{equation}

The detailed derivation from (2.52) to (2.58) can be found in Appendix B.

To evaluate the integral in (2.58), as shown in figure 2(b), we deform the integration path in (2.58) from $q/q_n$ to the point $\textrm {e}^{\pi \textrm {i}/2}$ , then along a unit circle to reach $q/q_n=1$ situated on the real axis, and finally along a ray from $q/q_n=1$ to $\infty$ . Using the identity (2.54), we can evaluate (2.58) to obtain

(2.59) \begin{equation} \begin{aligned} R_{+}(q)=\prod _{n=5}^8(q-q_n)^{1/2}\textrm {exp}\left(-\frac {\textrm i}{\pi }\sum _{n=1}^{8}\xi _n q_n \textrm {Li}_2(1-q/q_n) - \sum _{n=5}^{8}{\xi _n q_n}\textrm {ln}_+ \frac {q}{q_n}\right). \end{aligned} \end{equation}

Here, $\textrm {Li}_2$ represents the polylogarithm function of the second order (Gradshteyn & Ryzhik Reference Gradshteyn and Ryzhik1980).

Using (2.46), (2.48) and (2.59), we obtain the final form of $K_{+}$ , i.e.

(2.60) \begin{equation} K_{+}(k)=\frac {\sqrt {k-M_2}}{T_{+}(k)R_{+}(k/S^{3/4})\mathcal {Q}_{+}(k/S^{3/4})}. \end{equation}

A complete expression of $K_{+}$ is

(2.61) \begin{equation} \begin{aligned} K_{+}(k)=\sqrt {\pi }&\frac {\Gamma \left(1+\sqrt {M_{\textrm {i}}^2-1}k/2\pi \right)}{\Gamma \left(1/2+\sqrt {M_{\textrm {i}}^2-1}k/2\pi \right)}\frac {\sqrt {k-M_2}\sqrt {k/S^{3/4}+0}}{\prod _{n=5}^8(q-q_n)^{1/2}} \\ &\textrm {exp}\left(\frac {\textrm i}{\pi }\sum _{n=1}^{8}\xi _n q_n \textrm {Li}_2(1-q/q_n) + \sum _{n=5}^{8}{\xi _n q_n}\textrm {ln}_+ \frac {q}{q_n}\right). \end{aligned} \end{equation}

As $k\rightarrow 0$ , the Taylor expansion of $K_{+}(k)$ reduces to

(2.62) \begin{equation} K_{+}(k)\approx \frac {\textrm {i}\pi \sqrt {M_2}}{\mathcal {C}S^{3/8}}{\sqrt {k}}+O\left(k^{3/2}\right). \end{equation}

Here, the constant $\mathcal {C}$ is defined as $\mathcal {C}=\prod _{n=5}^{8}q_n^{1/2}$ . Equation (2.62) will be used to determine $\epsilon$ in § 2.3.

2.3. Instability waves excited by the acoustic waves at the nozzle lip

Under the excitation of the acoustic wave, new instability waves would emerge near the nozzle lip (Crighton Reference Crighton1985). We can determine this instability wave by imposing the unsteady Kutta condition at the nozzle lip (Crighton Reference Crighton1981). Following Crighton (Reference Crighton1972) and considering that we focus on the upper half-plane ( $y\gt 1/2$ ), we assume that the excited instability wave can be represented by

(2.63) \begin{equation} \phi ^{*}=\mathcal {A}\textrm {e}^{-\textrm {i}\mathcal {K}x-\gamma (\mathcal {K})y}, \end{equation}

where $\mathcal {A}$ represents the amplitude of the instability wave and $\gamma (\mathcal {K})$ is defined by $\gamma (\mathcal {K})=\sqrt {\mathcal {K}^{2}-\omega ^2 M_{\textrm {o}}^2}$ . The parameter $\mathcal {K}$ denotes the wavenumber of the instability waves of the antisymmetric mode and its determination is outlined in Appendix A. The choice of the form in (2.63) is believed to be reasonable, particularly considering the sinuous nature exhibited by both the upstream forcing and the scattered acoustic field.

The instability wave then propagates downstream and subsequently interacts with the shock structures (Powell Reference Powell1953b ; Suda et al. Reference Suda, Manning and Kaji1993; Edgington-Mitchell et al. Reference Edgington-Mitchell, Wang, Nogueira, Schmidt, Jaunet, Duke, Jordan and Towne2021). This interaction leads to the emission of sound, which serves as the upstream forcing and completes the screech feedback loop. To determine the instability amplitude $\mathcal {A}$ in (2.63), we first calculate the scattered field resulting from the newly excited instability wave at the nozzle lip. The singular component within the resulting scattered field needs to be eliminated using the singularity in (2.38). This process fixes the value of $\mathcal {A}$ to satisfy the unsteady Kutta condition.

Replacing the forcing term $GY$ in (2.9) and (2.10) with $\phi ^{*}$ while following the same procedures outlined in § 2.1, we obtain the scattered field due to the presence of the instability wave, i.e.

(2.64) \begin{equation} \phi ^{*^\prime }_{+}(k,1/2)=\frac {\textrm {i}\mathcal {A}\gamma (\mathcal {K})\textrm {e}^{-\gamma (\mathcal {K})/2}}{k-\mathcal {K}}\left (1-\frac {K_{+}(k)}{K_{+}(\mathcal {K})}\right ). \end{equation}

The component of $\phi ^{*^\prime }_{+}(k)$ similar to (2.38) that leads to singularity is

(2.65) \begin{equation} -\frac {\textrm {i}\mathcal {A}\gamma (\mathcal {K})\textrm {e}^{-\gamma (\mathcal {K})/2}}{K_{+}(\mathcal {K})}\left (\frac {K_{+}(k)}{k}\right ). \end{equation}

The unsteady Kutta condition demands no singularity in the vicinity of $r\rightarrow 0$ ; therefore, from (2.38) and (2.65), we must have

(2.66) \begin{equation} \frac {\textrm {i}G}{K_{+}(\textrm {i}\epsilon )} -\frac {\textrm {i}\mathcal {A}\gamma (\mathcal {K})\textrm {e}^{-\gamma (\mathcal {K})/2}}{K_{+}(\mathcal {K})}=0. \end{equation}

The parameter $G$ and the kernel function $K_{+}$ are defined by (2.7) and (2.60), respectively. To evaluate $\mathcal {A}$ , we need to determine the small parameter $\epsilon$ , which seems to be an arbitrary value introduced to address the issue of non-convergence. However, $\epsilon$ in fact has a physics-imposed value and we may uniquely determine it by re-evaluating the problem without approximations.

First, replace the forcing term $GY$ in (2.9) and (2.10) by

(2.67) \begin{equation} F(h)=\frac {1}{2\pi }\int _{-\infty }^{\infty }D_{{\textrm {o}}}(\lambda )\textrm {exp}(-\textrm {i}\lambda h)\textrm {d}\lambda . \end{equation}

Then the Wiener–Hopf equation (2.30) takes a new form, i.e.

(2.68) \begin{equation} \phi _{+}^{\prime }(k,1/2)+\mathcal {F}(k)=K(k)F_{-}(k). \end{equation}

Here,

(2.69) \begin{equation} \begin{aligned} \mathcal {F}(k) &= \int _{-\infty }^{0} \int _{0}^{\infty } \mathcal {G}(\zeta ^{*})\textrm {exp}(\zeta ^{*} x) \textrm {e}^{\textrm {i}kx} \textrm {d}\zeta ^{*} \textrm {d}x, \end{aligned} \end{equation}

where $\zeta ^{*}=-\textrm {i}\lambda$ , and the form of $\mathcal {G}(\zeta ^{*})$ is shown in Appendix B. If we intend to completely account for the effect of the upstream forcing, we can replace $G$ by $\mathcal {G}(\zeta ^{*})$ and $\epsilon$ by $\zeta ^{*}$ , and then carry out the integration from $0$ to $\infty$ . The mathematical proof of this procedure can be found in Appendix C. Thus, it follows that

(2.70) \begin{equation} \frac {G}{K_{+}(\textrm {i}\epsilon )}=\int _{0}^{\infty } \frac {\mathcal {G}(\zeta ^{*})}{K_{+}(\textrm {i}\zeta ^{*})} \textrm {d}\zeta ^{*}. \end{equation}

Figure 4.

Value of the parameter $\epsilon$ in the range of $1.2\lt M_{\textrm {i}}\lt 2$ at the screeching frequency. The frequency is obtained by the original formula (1.1) proposed by Powell (Reference Powell1953a ).

The integration on the right-hand side of (2.70) can be evaluated using Watson’s lemma. On the left-hand side of (2.70), given that $\epsilon$ is a small term, $K_{+}(\textrm {i}\epsilon )$ can be evaluated employing its Taylor expansion. By combining (2.62) and (C7), we have

(2.71) \begin{equation} \epsilon =-\frac {\textrm {i}G^2h^3}{\pi \mathcal {G}^2(0)}. \end{equation}

To verify that $\epsilon$ is a small term, we calculate it using (2.71) and plot it in figure 4. The angular frequency of the instability wave is calculated using the original formula (1.1) proposed by Powell (Reference Powell1953a ). Figure 4 shows that within the range of $1.2\lt M_{\textrm {i}}\lt 2$ , the value of $\epsilon$ initially decreases from its maximal value of $5.9\times 10^{-3}$ at $M_{\textrm {i}}=1.2$ to $4.9\times 10^{-4}$ at $M_{\textrm {i}}=1.6$ , and then gradually increases to approximately $1.9\times 10^{-3}$ at $M_{\textrm {i}}=2$ . In the entire Mach number region, $\epsilon$ is small, ensuring the validity of using Taylor expansion in (2.62) and the omission of $G\epsilon /\omega$ in (2.25).

Having obtained all the parameters, the instability wave is found to be

(2.72) \begin{equation} \phi ^{*}(0, 1/2)=\frac {G K_{+}(\mathcal {K})}{\gamma (\mathcal {K})K_{+}(\textrm {i}\epsilon )}. \end{equation}

This equation enables us to define the important transfer function, i.e.

(2.73) \begin{equation} {H}(\omega , M_{\textrm {i}}, \nu )=\frac {K_{+}(\mathcal {K})}{\gamma (\mathcal {K})K_{+}(\textrm {i}\epsilon )}, \end{equation}

where $K_{+}$ is given by (2.60), $\mathcal {K}$ by (A7), $\gamma (\mathcal {K})=\sqrt {\mathcal {K}^{2}-\omega ^2 M_{\textrm {o}}^2}$ and $\epsilon$ by (2.71). This transfer function describes the linear mapping between the newly excited instability waves and the upstream propagating acoustic waves. For clarity, Appendix D provides a summary of the established models upon which the present model is built, the models used for validation and the underlying assumptions for these models. Additionally, the fundamental assumptions employed in developing the transfer function are also outlined in Appendix D. Using this function, both the scattering efficiency and the phase delay between these two waves can be readily examined.