2.1 Base flow

The base flow behind a wing is a complex flow field for which no analytical solution exists; therefore, the base flow is determined numerically. We utilize direct numerical simulation to compute the flow around a NACA0012 profile with a flat wingtip positioned at an angle of attack of $5^{\circ }$ . The chord Reynolds number is fixed at $Re_{c}=1000$ with a half-span $b/2=1.25$ and a chord $c=1$ , allowing for all spatial coordinates to correspond to their non-dimensional values based on the chord length. In the transverse ( $y$ ) direction, the computational domain extends over the interval $y\in [-15,\,15]$ , while in the spanwise ( $z$ ) direction we take $z\in [-1.25,\,14.375]$ . The streamwise direction extends in a semi-circular arc with a radius of 15 and progresses downstream 15 chords, discretized by a C-grid.

An incompressible, finite-volume flow solver, Cliff (CharLES package), developed by Cascade Technologies (Ham & Iaccarino Reference Ham and Iaccarino2004; Ham, Mattsson & Iaccarino Reference Ham, Mattsson and Iaccarino2006), is used, which is second-order accurate in both time and space. The boundary conditions at the inlet are prescribed with $(U,V,W)=(1,0,0)$ , and at the outlet a convective outflow boundary condition is employed. In the far field and on the symmetry plane, we impose free-slip boundary conditions.

The flow around the wing is simulated with approximately 3.5 million elements. The spatial field is discretized on a hybrid mesh, where the boundary layer and wake regions are covered by a structured grid, while the far field outside this region is computed on an unstructured grid consisting of tetrahedral elements. To ensure the grid is sufficient, we locally increase the number of elements in the wake by approximately 2.1 million in the wake region, showing that the refined case recreated the same flow field that contains minimal differences. The computational domain is chosen such that the cross-section of the wing contains a blockage ratio of less than 0.5 %, indicating that the domain size is sufficiently large. These results provide us with confidence that the base flow accurately captures the trailing-line vortex in free space.

The simulation ran until it reached a steady-state base flow, shown in figure 2. The blue iso-surface visualizes the trailing vortex with the $Q$ -criterion (Jeong & Hussain Reference Jeong and Hussain1995). As shown, the vortex core evolves smoothly downstream, with the radius slightly diminishing in size with downstream progression due to viscous diffusion. A contour slice shows the streamwise velocity deficit in the wake of the airfoil at $x=3$ , i.e. the plane to be examined in this paper. Our study focuses on this plane because flow at this location is influenced by both the wake and the trailing vortex allowing us to gain insight into both wake and trailing vortex instabilities.

Figure 2. The base flow at $Re_{c}=1000$ , visualized with blue iso-surfaces, showing the $Q$ -criterion ( $Q=0.25$ ). The dark grey region corresponds to the wing location. The contour slice of the streamwise velocity at $x=3$ illustrates the streamwise velocity deficit in the wake region, while the $Q$ -criterion visualizes the structure of the trailing vortex.

The downstream development of the streamwise ( $U$ ), transverse ( $V$ ) and spanwise ( $W$ ) velocity profiles are shown in figure 3 at $z=-1.25$ and $-0.25$ (i.e. midspan of the wing and 80 % span, respectively). At the midspan of the wing (corresponding to $z=-1.25$ , figure 3 a), a streamwise velocity deficit is observed near the trailing edge. This deficit advects downward, indicated by the dash-dotted line of the minimum velocity, due to the downwash caused by the induced negative transverse velocity. At $z=-0.25$ (figure 3 b), the higher level of transverse velocity advects the streamwise velocity deficit farther upstream, causing the wake to descend more rapidly. The spanwise velocity, $W$ , is zero at the midspan location, but increases approximately to the level of the transverse velocity at $z=-0.25$ . At $x=3$ , the effect of the wake is significant, but not so overwhelming as to overshadow the effect of the developing vortex on the stability analysis. Hence, we select $x=3$ as a representative location to examine the combined efforts from the wake and trailing vortex onto the stability of the flow.

Figure 3. Streamwise ( $U$ , dashed line), transverse ( $V$ , solid line) and spanwise ( $W$ , dotted line) velocities ( $V$ and $W$ are scaled up by a factor of 5 for graphical clarity) of the wake of the wing at $z=-1.25$ (a) and $z=-0.25$ (b) with downstream progression from $x=0$ to 5. The dash-dotted line with the solid-dot markers shows the point of minimum axial velocity deficit in the wake. The $x$ -locations are $x=0,1.5,3$ and 4.5.

2.2 Stability analysis approach and validation

We examine the modal behaviour of small perturbations about the given base flow. Following standard linear stability analysis formulation, the state variable can be decomposed into the steady base flow and an unsteady disturbance as

(2.1) $$\begin{eqnarray}\displaystyle \tilde{\boldsymbol{q}}(x,y,z,t)=Q(x,y,z)+\unicode[STIX]{x1D716}\hat{q}(x,y,z,t). & & \displaystyle\end{eqnarray}$$

Here, $\tilde{\boldsymbol{q}}=(\tilde{u} ,\,\tilde{v},\,\tilde{w},\,\tilde{p})^{\text{T}}$ denotes the state vector, $Q=(U,\,V,\,W,\,P)^{\text{T}}$ is the known steady base flow, $\hat{q}=(\hat{u} ,\,\hat{v},\,{\hat{w}},\,\hat{p})^{\text{T}}$ represents the disturbance and $\unicode[STIX]{x1D716}$ stands for a small disturbance amplitude. By only considering small disturbances, we neglect the nonlinear, $O(\unicode[STIX]{x1D716}^{2})$ terms in the incompressible Navier–Stokes equations, which in turn results in the linearized disturbance equations (Schmid & Henningson Reference Schmid and Henningson2001). In the intermediate region behind the airfoil, the slow process of viscous diffusion relative to streamwise convection produces a quasi-parallel flow, with streamwise gradients being less than 2.75 % of the transverse gradients. This quasi-parallel assumption on streamwise variations of the base flow leads to the disturbance equations with constant coefficients in the streamwise direction and permits wave-like disturbance solutions of the form

(2.2) $$\begin{eqnarray}\displaystyle \hat{q}(x,y,z,t)=q(y,z)\exp (\text{i}\unicode[STIX]{x1D6FC}x-\text{i}\unicode[STIX]{x1D714}t). & & \displaystyle\end{eqnarray}$$

In the above expression, $q(y,z)$ denotes the shape function, and the argument in the exponent of (2.2) is the phase function, where $\unicode[STIX]{x1D6FC}$ and $\unicode[STIX]{x1D714}$ represent the streamwise wavenumber and radian frequency, respectively. Substitution of (2.2) into the linearized disturbance equations results in a generalized eigenvalue problem. The specific form of the eigenvalue problem depends on the type of stability analysis. For a temporal stability analysis, we assume growth/decay in time and allow for a complex frequency $\unicode[STIX]{x1D714}=\unicode[STIX]{x1D714}_{r}+\text{i}\unicode[STIX]{x1D714}_{i}$ , with $\unicode[STIX]{x1D714}_{i}>0$ denoting exponential growth at the radian frequency $\unicode[STIX]{x1D714}_{r}$ . The associated eigenvalue problem reads ${\mathcal{A}}_{T}\boldsymbol{q}=\unicode[STIX]{x1D714}{\mathcal{B}}_{T}\boldsymbol{q}$ where ${\mathcal{A}}_{T}$ and ${\mathcal{B}}_{T}$ are defined in appendix A. Conversely, for a spatial stability analysis, the streamwise wavenumber $\unicode[STIX]{x1D6FC}$ represents the complex eigenvalue $\unicode[STIX]{x1D6FC}=\unicode[STIX]{x1D6FC}_{r}+\text{i}\unicode[STIX]{x1D6FC}_{i}$ with $\unicode[STIX]{x1D6FC}_{i}<0$ indicating spatial exponential growth. The corresponding eigenvalue problem is nonlinear in the eigenvalue $\unicode[STIX]{x1D6FC},$ but can be recast into an equivalent linear problem ${\mathcal{A}}_{S}\boldsymbol{q}_{S}=\unicode[STIX]{x1D6FC}{\mathcal{B}}_{S}\boldsymbol{q}_{S}$ via a companion matrix technique (Tisseur & Meerbergen Reference Tisseur and Meerbergen2001). The operators ${\mathcal{A}}_{S}$ and ${\mathcal{B}}_{S}$ are also defined in appendix A.

The spatial discretization of the respective stability matrices ${\mathcal{A}}$ and ${\mathcal{B}}$ is performed using the Chebyshev-based spectral collocation method (Kopriva Reference Kopriva2009). The Chebyshev Gauss–Lobatto points $\unicode[STIX]{x1D702}_{j}$ are mapped onto the transverse coordinate direction $y$ using (Hein & Theofilis Reference Hein and Theofilis2004)

(2.3) $$\begin{eqnarray}\displaystyle y_{j}=g(\unicode[STIX]{x1D702}_{j};\unicode[STIX]{x1D705},y_{\infty })=y_{\infty }\frac{\tan \left(\displaystyle \frac{\unicode[STIX]{x1D705}\unicode[STIX]{x03C0}\unicode[STIX]{x1D702}_{j}}{2}\right)}{\tan \left(\displaystyle \frac{\unicode[STIX]{x1D705}\unicode[STIX]{x03C0}}{2}\right)},\quad j=0,1,\ldots ,N_{y}, & & \displaystyle\end{eqnarray}$$

where $y_{\infty }$ denotes the largest far-field $y$ -value noted below, and the parameter $\unicode[STIX]{x1D705}$ adjusts the spread of the points about $y=0$ . When examining the discrete branches that are isolated in the vortex regions, we use a value of $\unicode[STIX]{x1D705}=0.96$ . In contrast, when analysing the continuous branch, we require a larger number of points in the free stream, and thus relax the parameter to $\unicode[STIX]{x1D705}=0.85$ , yielding more gradual stretching of the grid.

In the spanwise coordinate direction with $z\in [-1.25,14.375]$ , two nested mappings are utilized. First, equation (2.3) is used to map $\unicode[STIX]{x1D701}_{j}\rightarrow \unicode[STIX]{x1D701}_{j}^{\ast }$ , where $\unicode[STIX]{x1D701}_{j}$ denotes the Gauss–Lobatto points in the $z$ -direction, and $\unicode[STIX]{x1D701}_{j}^{\ast }$ are the mapped points per (2.3). Throughout this study, $\unicode[STIX]{x1D705}$ in the $z$ -direction is set to be equal to the value used in the $y$ -direction. We maintain the far-field value to be $\unicode[STIX]{x1D701}_{\infty }^{\ast }=1$ to preserve the scale on $\unicode[STIX]{x1D701}^{\ast }\in [-1,1]$ . With the points reclustered about the centre, we then apply the second mapping from $\unicode[STIX]{x1D701}^{\ast }$ to $z\in [-1.25,14.375]$ using a rational function (Hanifi, Schmid & Henningson Reference Hanifi, Schmid and Henningson1996)

(2.4) $$\begin{eqnarray}\displaystyle z_{j}=h(\unicode[STIX]{x1D701}_{j}^{\ast };z_{min},z_{mid},z_{\infty })=z_{min}+a\frac{1+\unicode[STIX]{x1D701}_{j}^{\ast }}{b-\unicode[STIX]{x1D701}_{j}^{\ast }},\quad j=0,1,\ldots ,N_{z}, & & \displaystyle\end{eqnarray}$$

with

(2.5a,b ) $$\begin{eqnarray}\displaystyle a=\frac{z_{mid}z_{\infty }}{z_{\infty }-2z_{mid}},\quad b=1+\frac{2a}{z_{\infty }}. & & \displaystyle\end{eqnarray}$$

We choose $z_{min}=-1.25$ , $z_{mid}=1.25$ and $z_{\infty }=15.625$ . The above expression distributes half the mapped collocation points to cluster below and above $z_{mid}$ prior to translating the data by $z_{min}$ .

Stability analyses show relative insensitivity of the leading eigenvalue to the choice of domain size owing to the exponential decay of its eigenfunction outside the shear layer. For this reason, final computations are performed on the domain of $y\in [-10,10]$ and $z\in [-1.25,10]$ in order to achieve satisfactory resolution within the vortex and wake regions. The discretized eigenvalue problem is solved using the implicitly restarted Arnoldi method with the shift-and-invert technique implemented in ARPACK (Lehoucq, Sorensen & Yang Reference Lehoucq, Sorensen and Yang1996).

Table 1. Validation of the present temporal and spatial analyses with duct flow (Theofilis, Duck & Owen Reference Theofilis, Duck and Owen2004) and the Batchelor vortex (Paredes Reference Paredes2014), respectively. Shown are the least stable eigenvalues with convergence. Also included (rightmost column) are results of our temporal stability analysis of the flow past a wingtip for $\unicode[STIX]{x1D6FC}=5.5$ at $Re_{c}=1000$ , demonstrating convergence of the converged least stable eigenvalue as the grid is refined.

Our bi-global stability code is validated against stability analyses of two-dimensional channel flow (Theofilis et al. Reference Theofilis, Duck and Owen2004; Paredes Reference Paredes2014) and the Batchelor vortex (Fabre & Jacquin Reference Fabre and Jacquin2004; Paredes Reference Paredes2014). The results of these validations are displayed in table 1. For the present analysis, convergence towards a relative error of $10^{-4}$ is achieved for a resolution of $N_{z}\times N_{y}=60\times 60$ . Therefore, for cases that require a large number of eigenvalue problems to be solved, we employ $N_{y}\times N_{z}=60\times 60$ , and for single case studies, we increase the number to $N_{y}\times N_{z}=80\times 80$ for more accurate results.

3.1 Parametric dependence on the streamwise wavenumber

Starting with a base-flow profile extracted at the streamwise position of $x=3$ , we solve the bi-global temporal stability problem and vary the streamwise wavenumber from $\unicode[STIX]{x1D6FC}=2$ to 6.5 in increments of $\unicode[STIX]{x0394}\unicode[STIX]{x1D6FC}=0.5$ . Recall that the $x=3$ location was chosen to examine both wake and tip vortex influences. The results are shown in figure 4. Three particular eigenvalues, representing a specific type of instabilities (see more details in § 3.3), are tracked: the first is termed the ‘principal wake instability’, the second eigenvalue is linked to a ‘vortex instability’ and the third, labelled ‘wake instability’, displays features consistent with a higher-order wake mode. The principal wake instability remains the primary instability throughout the range of $\unicode[STIX]{x1D6FC}$ examined, while the relative strength of the wake and vortex instabilities are dependent on the streamwise wavenumber. For low-wavenumber disturbances (i.e. large wavelengths), the vortex instability grows more rapidly than the wake instability.

Figure 4. Temporal growth rate, $\unicode[STIX]{x1D714}_{i}$ , at $x=3$ with varying streamwise wavenumber $\unicode[STIX]{x1D6FC}$ , with $N_{z}=N_{y}=60$ . The three instabilities are labelled as the principal wake instability, wake instability and vortex instability, according to their modal features. The dashed line indicates the wavenumber that we examine throughout the temporal analysis.

The maximum growth rate of the principal wake instability is $\unicode[STIX]{x1D714}_{i}\approx 0.34$ , which is achieved at a streamwise wavenumber of $\unicode[STIX]{x1D6FC}=5.5$ . At this wavenumber, both the wake and vortex instabilities show a growth rate of $\unicode[STIX]{x1D714}_{i}\approx 0.18$ and $\unicode[STIX]{x1D714}_{i}\approx 0.14$ , respectively, suggesting a balance between the vortex instability and the wake instability. With the wavenumber fixed at $\unicode[STIX]{x1D6FC}=5.5$ , we now present the bi-global spectrum and analyse the associated modes in detail.

3.2 The temporal bi-global spectrum – an overview

For a streamwise wavenumber of $\unicode[STIX]{x1D6FC}=5.5$ and a chord-based Reynolds number of $Re_{c}=1000$ , the temporal eigenvalue spectrum is displayed in the complex $\unicode[STIX]{x1D714}$ -plane in figure 5. The unstable half-plane $\unicode[STIX]{x1D714}_{i}>0$ is shaded in grey; instabilities, if any, lie in this upper half-plane. The quantitative values of these growth rates, $\unicode[STIX]{x1D714}_{i}$ , are approximate due to the steady-state nature of the base flow, implying the discrete branches may shift upward or downward in figure 5 (Sipp & Lebedev Reference Sipp and Lebedev2007). Despite this, however, we retain confidence in the relative levels and qualitative structure of these modes. Dividing $\unicode[STIX]{x1D714}_{r}$ by the wavenumber $\unicode[STIX]{x1D6FC}$ yields the streamwise phase speed, i.e. the speed at which the associated modes travel. Invoking a critical-layer argument (Schmid & Henningson Reference Schmid and Henningson2001), this speed generally relates directly to the local speed of the base flow, which implies that the disturbances with $\unicode[STIX]{x1D714}_{r}<\unicode[STIX]{x1D6FC}$ are generally confined to regions of velocity deficit.

The bi-global spectrum shows distinct branches of discrete eigenvalues. The associated eigenmodes are spatially confined by the shear of the base flow which acts as a wave guide, containing the modes between the critical layers while causing exponential decay in the free stream. The discrete branch of the bi-global spectrum left of the continuous branch (magenta vertical line) contains two pertinent sub-branches, which we term a ‘wake branch’ and a ‘vortex branch’ (see figure 5). The disturbances of the wake branch are spatially localized in both the wake and vortex regions, usually along the entire span of the wake. The modal structures from the vortex branch are almost entirely restricted to the vortex region. The top of the discrete branch, containing the most unstable mode, represents a spatial structure that is predominantly localized in the wake, with minor contributions in the vortex core. We refer to this mode as a ‘principal wake instability’ (see figure 4). Progressing farther down the wake sub-branch towards larger decay rates (i.e. $\unicode[STIX]{x1D714}_{i}$ is decreasing), we observe increasing interaction between the wake and tip region, while the wake component shows structures with a higher spanwise wavenumber (see more details below). Furthermore, a discrete branch to the right of the continuous branch, termed the ‘azimuthal branch’, contains higher-order azimuthal modes ( $m>1$ ) localized solely in the vortex region and is discussed at the end of § 4.3. As these modes are stable, we reserve discussion of these modes until the spatial stability analysis in § 4. In addition to the discrete branch, we also detect a continuous branch comprising of free-stream oscillations that exponentially decay as they progress towards the wake and vortex of the base flow. The apex of the continuous branch lies at $\unicode[STIX]{x1D714}=\unicode[STIX]{x1D6FC}-\text{i}\unicode[STIX]{x1D6FC}^{2}/Re_{c}$ . Further details of the continuous spectrum are presented in § 3.4.

Figure 5. The eigenvalue spectrum from a temporal stability analysis at $x=3$ for $Re_{c}=1000$ , $\unicode[STIX]{x1D6FC}=5.5$ and $N_{z}=N_{y}=80$ . The grey region represents the unstable half-plane. The most unstable mode at $\unicode[STIX]{x1D714}=4.862+0.340\text{i}$ corresponds to a wake-dominated mode denoted the principal wake instability. The spectrum is annotated to help provide a physical description of each branch of the spectrum.

3.3 The discrete branches

The discrete branch of the temporal spectrum shows an exponential instability whose structure is shown in figure 6 and is termed the principal wake instability. Figure 6(a) indicates (by a red circle) the location of the eigenvalue within the bi-global spectrum, (b) shows iso-surfaces of the $Q$ -criterion of the perturbation velocity for $Q=\pm 0.01$ , visualized in blue and red, respectively. The contour slice of the streamwise $U$ -velocity of the base flow is projected at $x=1$ to provide reference to the location of the wake and vortex region. Figure 6(c) shows the streamwise vorticity of the eigenfunction, with the solid and dashed lines visualizing positive and negative vorticity, respectively.

Figure 6. Bi-global principal wake instability mode. (a) The eigenvalue spectrum plotted with the corresponding eigenvalue circled in red. (b) Shows iso-surfaces of the $Q$ -criterion of the perturbation velocity with $Q=\pm 0.01$ (blue and red, respectively). The contour slice orients the figure to indicate where the wake is located. (c) Shows the streamwise vorticity with solid line contours representing positive vorticity and dashed-lined contours denoting negative vorticity. For a movie of (c), see movie 1 in the supplemental material.

As the flow leaves the trailing edge, the wake of the wing contains a streamwise velocity deficit. This deficit confines the wake instabilities and reduces the phase speed of the disturbances in this region, with the velocity gradients of the base flow shearing the disturbance, as shown by the aft-angle defined by the red and blue iso-surfaces. Little spanwise variation is present in this mode until the vortex is reached. The rapid decay of the mode as the disturbance approaches the vortex suggests an inhibitive effect of the rotational motion on the wake mode; however, as the instability grows, a helical vortex component that co-rotates with the base flow (see e.g. movie 1 of the supplemental material available at https://doi.org/10.1017/jfm.2017.866) becomes apparent that establishes a coupling between the wake and the vortex. As shown from the figure 6(c), the perturbation is clearly localized to the base-flow shear, decaying exponentially outwards toward the free stream. Similar to many previously studied shear flows, the regions of shear act similar to walls, containing the mode within its critical layers and showing evanescent decay outside.

Progressing farther down the wake branch we encounter a second instability with $\unicode[STIX]{x1D714}=4.813+0.1908\text{i}$ , labelled the higher-order wake instability. The corresponding bi-global eigenfunction, shown in figure 7, contains disturbances localized in both regions of the wake and the vortex. The structure of this instability is similar to the principal wake instability (figure 6); however, it contains a spanwise structure with a higher wavenumber in this coordinate direction. The appearance of smaller scales in the spanwise direction gives rise to larger viscous diffusion and, consequently, a lower growth rate relative to the principal wake instability. Similar to the principal wake instability, the disturbance is localized in the region of shear and exponentially decays in the free stream (see figure 7 c). The vortex region just inboard of $(y,z)=(0,0)$ also shows a helical vortex structure that co-rotates with the base flow (see movie 2 of the supplemental material).

Figure 7. Bi-global wake instability mode. Variables are plotted following the same format as in figure 6. For a movie of (c), see movie 2 in the supplemental material.

Although the vortex instability, shown in figure 8, has the lowest growth rate ( $\unicode[STIX]{x1D714}_{i}=0.146$ ) of all unstable bi-global modes, its growth rate is comparable to the higher-order wake instability. This type of instability is localized in the vortex core, where the base-flow velocity is larger than in the wake, as reflected in a phase speed higher than for the other instabilities. Due to its localization in the vortex core, we hypothesize that the tip vortex drives this mode.

The vortex instability shows a two-lobed structure that co-rotates with the base flow. This pattern may be analogous to the helical mode of vortex instabilities. However, the azimuthal inhomogeneity of the base flow precludes an exact decomposition of the helical nature into Fourier modes (i.e. with azimuthal wavenumbers $m$ of integer values). Nonetheless, the observed co-rotation with the base flow hints at an azimuthal mode with $m=1$ (see movie 3 of the supplemental material), contrary to the Batchelor vortex which shows an instability for the $m=-1$ mode (Khorrami Reference Khorrami1991; Mayer & Powell Reference Mayer and Powell1992). This implies that the presence of the wake significantly alters the instability properties and furthermore necessitates one to cautiously treat the implication of an analysis of the trailing-line vortex in isolation from the wake. Although the disturbance is localized in the vortex region, the disturbance exists inboard of the vortex region, showing a slight coupling with the wake.

Figure 8. Bi-global vortex instability mode. Variables are plotted following the same format as in figure 6. For a movie of (c), see movie 3 in the supplemental material.

Continuing down the vortex branch, a stable vortex mode is shown in figure 9. This mode again shows a two-lobed structure that co-rotates with the base flow, but also contains higher-order structures inboard and beneath the vortex region, reminiscent of the trend along the wake branch. As time increases, the disturbance convects inboard from the tip region with a co-rotating helical mode, as shown in movie 4 of the supplemental material. For eigenfunctions farther down the vortex branch, not shown for brevity, higher-order structures appear in the vortex and a more pronounced interaction with the wake develops, together with increasingly smaller spanwise modal scales.

Figure 9. Bi-global stable vortex mode. Variables are plotted following the same format as in figure 6. For a movie of (c), see movie 4 in the supplemental material.

Farther along the wake branch, shown in figure 10, the slow-phase-speed mode remains restricted to the wake region. A spanwise variation of the mode is clearly visible, showing three lobes across the span. The modal component in the vortex region again represents a two-lobed structure akin to a helical mode that co-rotates with the base flow (see movie 5 of the supplemental material). Progressing farther down the branch (not shown), we encounter an increase in spanwise structures, while the helical nature of the vortex region remains virtually constant.

Figure 10. Bi-global stable wake mode. Variables are plotted following the same format as in figure 6. For a movie of (c), see movie 5 in the supplemental material.

From the analysis of the discrete branches of the bi-global spectrum, some general observations are worth noting. The unstable region of the discrete branch decomposes into two sub-branches which capture (i) a coupled dynamics of the wake–vortex system, and (ii) the instabilities of the vortex core. The principal instability arises on the wake branch for a structure that is predominately wake dominated with a coupled component in the vortex region. The remaining two instabilities, from each of the branches, show similar growth rates, but differ in phase velocity. The vortical components of higher modes show two-lobed structures that co-rotate with the base-flow vorticity – in contrast to the stability behaviour of a Batchelor vortex that counter-rotates with the base flow. This discrepancy points towards a significant effect of the wake component on the dynamics of the trailing vortex in this intermediate region. We hypothesize that the presence of the wake aft of the wing imparts preference to the vortex modes to co-rotate with the base flow, shown in movies 1–5 of the supplemental material. Advancing along both branches towards larger decay rates, an increase in spanwise structures is observed in the wake, as is a more pronounced coupling between the wake and vortex dynamics.

3.4 The continuous branch

While the discrete part of the bi-global spectrum is dictated by the presence of shear in the base-flow velocity field, the fact that we consider viscous flow in an infinite domain introduces a continuous spectrum. This phenomenon has long been recognized and studied in simpler flows such as the Blasius boundary layer (Mack Reference Mack1976; Grosch & Salwen Reference Grosch and Salwen1978). In order to study the continuous branch, we relax the boundary condition at infinity and allow disturbances that are merely bounded at the computational boundary, and thus oscillate in the free stream. Physically, this part of the bi-global spectrum contains information about the interaction characteristics of free-stream perturbations with their discrete counterparts, and thus describes receptivity processes (Saric, Reed & Kerschen Reference Saric, Reed and Kerschen2002). The dispersion relation associated with the continuous spectrum gives insight into the filter behaviour of the flow to external perturbations. This type of analysis uncovers which far-field perturbations interact with the discrete instabilities and which perturbations do not pass the shear region to trigger modal growth.

Following standard procedure (see e.g. Schmid & Henningson Reference Schmid and Henningson2001), we derive an analytical expression for the continuous spectrum by taking the limit as $|y|,|z|\rightarrow \infty$ , eliminating base-flow gradients and allowing the approximation of $(U,V,W)=(U_{\infty },0,0)$ . This step reduces the governing equations to a system of partial differential equations with constant coefficients and thus allows solutions of the form $\boldsymbol{v}=\hat{\boldsymbol{v}}\exp (\text{i}\unicode[STIX]{x1D6FD}y+\text{i}\unicode[STIX]{x1D6FE}z)$ with $\unicode[STIX]{x1D6FD}$ and $\unicode[STIX]{x1D6FE}$ denoting wavenumbers in the respective coordinate directions. Upon substitution, we obtain a system of algebraic equations for $\unicode[STIX]{x1D6FD}$ and $\unicode[STIX]{x1D6FE}$ . To further reduce this system of four equations, we eliminate the streamwise disturbance velocity and pressure by utilizing the continuity equation and streamwise momentum equation, while maintaining the disturbance velocities, $v$ and $w$ . This results in an eigenvalue problem of the form

(3.1) $$\begin{eqnarray}\displaystyle \hat{{\mathcal{A}}}\hat{\boldsymbol{v}}=\unicode[STIX]{x1D714}\hat{{\mathcal{B}}}\hat{\boldsymbol{v}}, & & \displaystyle\end{eqnarray}$$

where $\hat{{\mathcal{A}}}$ and $\hat{{\mathcal{B}}}$ are derived in appendix B, and $\hat{\boldsymbol{v}}=(\hat{v},{\hat{w}})^{\text{T}}$ . The eigenvalues then are determined as solutions of a quadratic equation according to

(3.2) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D714}=\frac{-B\pm \sqrt{B^{2}-4AC}}{2A}, & & \displaystyle\end{eqnarray}$$

with

(3.3a-c ) $$\begin{eqnarray}\displaystyle A=\unicode[STIX]{x1D6FC}^{2}(\unicode[STIX]{x1D709}^{2}-\unicode[STIX]{x1D6FC}^{2}),\quad B=2\text{i}{\hat{c}}_{\ast }A,\quad C={\hat{c}}_{\ast }^{2}\unicode[STIX]{x1D6FC}^{2}(\unicode[STIX]{x1D6FC}^{2}-\unicode[STIX]{x1D709}^{2}), & & \displaystyle\end{eqnarray}$$

where

(3.4) $$\begin{eqnarray}\displaystyle {\hat{c}}_{\ast }=\text{i}\unicode[STIX]{x1D6FC}U_{\infty }+\frac{1}{Re_{c}}(\unicode[STIX]{x1D6FC}^{2}+\unicode[STIX]{x1D709}^{2}),\quad \unicode[STIX]{x1D709}^{2}=\unicode[STIX]{x1D6FD}^{2}+\unicode[STIX]{x1D6FE}^{2}. & & \displaystyle\end{eqnarray}$$

Restricting ourselves to real wavenumbers, i.e. $\unicode[STIX]{x1D6FD},\unicode[STIX]{x1D6FE}\in \mathbb{R}$ , yields an expression for the classical continuous branch: a line that extends from $\unicode[STIX]{x1D714}=\unicode[STIX]{x1D6FC}-(\unicode[STIX]{x1D6FC}^{2}/Re)\,\text{i}$ to $\unicode[STIX]{x1D714}=\unicode[STIX]{x1D6FC}-\infty \,\text{i}$ . This solution of the continuous branch is shown as a red line in figure 11. The fact that the continuous spectrum is a line, despite a dependence on the two parameters $\unicode[STIX]{x1D6FD}$ and $\unicode[STIX]{x1D6FE}$ , comes from the observation that $\unicode[STIX]{x1D714}$ depends on $\unicode[STIX]{x1D709}=\sqrt{\unicode[STIX]{x1D6FD}^{2}+\unicode[STIX]{x1D6FE}^{2}}$ not on $\unicode[STIX]{x1D6FD}$ and $\unicode[STIX]{x1D6FE}$ individually and thus is only parametrized by a single (composite) wavenumber.

Figure 11. Temporal continuous branch for $Re_{c}=1000$ and $\unicode[STIX]{x1D6FC}=5.5$ for $N_{y}=N_{z}=80$ . The red line corresponds to the theoretical continuous branch, while the symbols $\cdot$ correspond to the numerically computed eigenvalues from the bi-global stability analysis.

The discrete representation of the continuous branch of the bi-global spectrum appears to cover a two-dimensional wedge-like area. Modes that deviate from the analytical continuous branch show exponential decay as they approach the far-field boundary, an observation that has previously been reported (see, e.g. Obrist & Schmid Reference Obrist and Schmid2003a ,Reference Obrist and Schmid b ; Mao & Sherwin Reference Mao and Sherwin2011). Through inspection of the eigenfunctions, we conclude that these modes may be modelled as wavepackets (see § 3.5 for more details).

Eigenfunctions corresponding to eigenvalues along the continuous spectrum have support in the free stream and show oscillations until they reach the computational boundary. Moving down the continuous branch increases the number of oscillations in the free stream, as can also be deduced from the analytical expression (3.2). Modes farther from the continuous branch also oscillate in the free stream, see figure 12 for the eigenvalue $\unicode[STIX]{x1D714}=5.154-0.534\text{i}$ , but display an envelope that decays as $(y,z)\rightarrow (\pm \infty ,\infty )$ . As the modes deviate farther from the analytical continuous branch, the decay of the wavepacket envelope becomes increasingly rapid.

Figure 12. Eigenmode associated with the eigenvalue $\unicode[STIX]{x1D714}=5.154-0.534\text{i}$ , which lies off the continuous branch. As in figure 6 with $Q=0.005$ . With the eigenvalue off the continuous branch, the disturbance decays exponentially as it approaches the computational boundary, taking the form of a wavepacket.

To further explore this phenomenon, we relax the boundedness condition, used in the derivation of the continuous spectrum, to allow for exponential decay in the free stream. This is equivalent to permitting complex transverse wavenumbers (i.e. $\unicode[STIX]{x1D6FD},\unicode[STIX]{x1D6FE}\in \mathbb{C}$ ), where the real component corresponds, as before, to the free-stream oscillations, while the imaginary part describes the growth or decay in the transverse directions. This formulation allows for a rudimentary model of wavepackets. Since the disturbances must be bounded, we take the sign of the imaginary part appropriately to enforce exponential decay as we approach the free stream, $(y,z)\rightarrow (\pm \infty ,\infty )$ .

With the wavenumbers $\unicode[STIX]{x1D6FD}$ and $\unicode[STIX]{x1D6FE}$ allowed to take on complex values, we can map the lower half of the complex $(\unicode[STIX]{x1D6FD},\unicode[STIX]{x1D6FE})$ -plane under the analytical expression (3.2) for the continuous spectrum. Figure 13(b) shows the resulting parabolic spread in the complex $c$ -plane. With increasing magnitudes of $\unicode[STIX]{x1D6FD}_{i}$ and $\unicode[STIX]{x1D6FE}_{i}$ these parabolas spread upward and outward, and the mapped half-plane covers an increasingly larger area of the complex $c$ -plane. It is important to note that these parabolas cover a continuous area; the exact locations of the global eigenvalues in this region, however, are affected by the numerical discretization and the choice of computational domain size and depend sensitively on these numerical parameters.

Figure 13. The complex $(\unicode[STIX]{x1D6FD},\unicode[STIX]{x1D6FE})$ -plane shown in (a) is then mapped to the complex phase-speed plane, where $c=\unicode[STIX]{x1D714}/\unicode[STIX]{x1D6FC}$ shown in (b). The complex $(\unicode[STIX]{x1D6FD},\unicode[STIX]{x1D6FE})$ results in a continuous spectrum of an area in the complex $c$ -plane, rather than a line ( $\unicode[STIX]{x1D6FD}_{i},\unicode[STIX]{x1D6FE}_{i}=0$ ) as in classical stability theory.

3.5 Wavepacket analysis

We proceed by introducing a third manner of analysis of the spectral problem of trailing-line vortices. Our first analysis in § 3.3 solved the full global eigenvalue problem and was particularly suited to regions where the coefficients of our system (given by the base flow and its derivatives) are rapidly varying; this analysis resulted in the discrete modal structures. The second approach in § 3.4 considered the limit of large distances from the regions of shear, where the base flow appears uniform and our system can be approximated by a constant-coefficient set of equations; the resulting wave solutions and the corresponding analytic dispersion relation (3.2) constitute the continuous branch of the spectrum and describe oscillatory modal solutions that spatially continue to infinity. In this section, we address the intermediate regime where eigensolutions can be approximated by compact wavepackets. This approximation becomes exponentially accurate as a small, user-defined parameter, in our case $h=1/\sqrt{Re_{c}}$ , tends to zero. This type of analysis covers structures in the outer neighbourhoods of the vortex and wake and sheds light on the remaining global eigenvalues which can neither be identified as discrete nor as continuous. An analysis of this type has been successfully performed in simpler configurations (Obrist & Schmid Reference Obrist and Schmid2010; Mao & Sherwin Reference Mao and Sherwin2011); here, we extend it to a two-dimensional base flow. For the sake of focus, the mathematical details of the methodology are relegated to the appendix C, while the motivation, a conceptual summary, the results and physical interpretation are provided in this section.

Wavepacket analysis helps depict the two-dimensionality of the eigenvalue spectrum of the bi-global analysis through the solution of a system of algebraic equations. This methodology may yield efficient means of computing the spectrum that can be applied to application of transient growth and potential applications such as aeroacoustics where wavepackets are of prime importance. This motivates the extension of wavepacket analysis from the simplified one-dimensional cases shown in Trefethen (Reference Trefethen2005), Obrist & Schmid (Reference Obrist and Schmid2010), Mao & Sherwin (Reference Mao and Sherwin2011) to the two-dimensional case where more complex flows are applicable.

Conceptually, wavepacket analysis assumes eigenfunctions of the general form $\exp (\text{i}f(y,z;h))$ with $f$ denoting an unknown phase function that depends on both spatial coordinates as well as a small parameter $h$ . This approach is reminiscent and akin to a Wentzel–Kramers–Brillouin–Jeffreys (WKBJ) procedure and utilizing wavepacket analysis provides an approximate dispersion relation. Upon substitution of our general form into the governing equations and introducing the spatially local wavenumbers $\unicode[STIX]{x1D735}f=(f_{y},f_{z})^{\text{T}}\equiv (1/h)(\unicode[STIX]{x1D6FD},\unicode[STIX]{x1D6FE})^{\text{T}}$ , we arrive at a dispersion relation for $\unicode[STIX]{x1D714}$ that depends on both $(y,z)$ and the associated wavenumbers $(\unicode[STIX]{x1D6FD},\unicode[STIX]{x1D6FE})$ and thus contains a hybrid of physical and spectral components. This dispersion relation is approximate but approaches the true dispersion relation exponentially as $h\rightarrow 0$ . Mathematically, this approximate dispersion relation is referred to as the symbol (Trefethen Reference Trefethen2005). Thus far, solutions of the form $\exp (\text{i}f(y,z;h))$ do not necessarily represent wavepackets. Additional constraints have to be imposed to ensure compactness of the eigenfunctions. A simple second-order Taylor expansion about the peak of the wavepacket $\boldsymbol{x}_{\ast }=(y_{\ast },z_{\ast })$ yields (we will use the subscript $_{\ast }$ to denote evaluation at peak location)

(3.5) $$\begin{eqnarray}\displaystyle \exp (\text{i}f(y,z;h))\sim \underbrace{\exp (\text{i}f_{\ast })}_{const.}\times \underbrace{\exp (\text{i}\boldsymbol{b}_{\ast }^{\text{T}}(\boldsymbol{x}-\boldsymbol{x}_{\ast }))}_{carrier~wave}\times \underbrace{\exp (\text{i}(\boldsymbol{x}-\boldsymbol{x}_{\ast })^{\text{T}}\mathsf{C}_{\ast }(\boldsymbol{x}-\boldsymbol{x}_{\ast }))}_{packet~envelope}, & & \displaystyle\end{eqnarray}$$

with $\boldsymbol{b}=\unicode[STIX]{x1D735}f$ and $\mathsf{C}=\unicode[STIX]{x1D6FB}^{2}f$ as the gradient and Hessian of $f$ , respectively. Clearly, a condition on $\mathsf{C}$ has to be imposed to ensure that the wavepacket envelope is compact, rather than divergent. This condition is referred to as the twist condition (Trefethen Reference Trefethen2005) and is given in appendix C.

The procedural steps of a wavepacket analysis are then as follows. We sweep over $y_{\ast }$ and $z_{\ast }$ (the peak location of our wavepacket mode), covering the physical domain of interest, e.g. the neighbourhood of the vortex or the wake; simultaneously, we also sweep over the wavenumbers $\unicode[STIX]{x1D6FD}_{\ast }$ and $\unicode[STIX]{x1D6FE}_{\ast }$ , covering the spatial scales of interest. For each point $(y_{\ast },z_{\ast },\unicode[STIX]{x1D6FD}_{\ast },\unicode[STIX]{x1D6FE}_{\ast })$ we verify the twist condition and discard points that do not satisfy it, deeming them unphysical. The accepted points are then mapped, via the symbol (i.e. approximate dispersion relation), into the complex domain where they trace out a parameterized curve as $(y_{\ast },z_{\ast },\unicode[STIX]{x1D6FD}_{\ast },\unicode[STIX]{x1D6FE}_{\ast })$ are varied. Each point on this curve constitutes an (exponentially) approximate eigenvalue associated with our wavepacket eigenmode. Depending on the complexity of the symbol, these curves can intersect one another on the complex plane, giving rise to multiple wavepacket solutions. The parameterization by four parameters $(y_{\ast },z_{\ast },\unicode[STIX]{x1D6FD}_{\ast },\unicode[STIX]{x1D6FE}_{\ast })$ generally yields an area in the complex plane that constitutes the spectral location of approximate wavepacket solutions to our stability problem. The wedge-like clustering of discrete points about the continuous spectrum in figure 5 about $\unicode[STIX]{x1D714}_{r}\approx 5.5$ is a (discrete) manifestation of this area spectrum. For more details, refer to appendix C.

The symbol can be used to assess the characteristics of the bi-global eigenvalue spectrum. In figure 14(a), we select an eigenvalue (circled in red) to compare the true eigenmode from our bi-global analysis with the corresponding approximate mode computed through wavepacket analysis. To this end, we determine the parameters  $y_{\ast }$ and  $z_{\ast }$ , i.e. the location of the wavepacket, using (C 5) and $\unicode[STIX]{x1D6FD}_{\ast }/h$ , $\unicode[STIX]{x1D6FE}_{\ast }/h,$ i.e. the approximate wavenumbers, using (C 4). We then substitute these parameters into (C 5) to obtain four distinct wavepackets. Although this provides the location and shape of the wavepackets, the weighting of each individual approximation is still unknown. We use the eigenfunction from the bi-global spectrum and minimize the least-squares error between a linear combination of approximations and the true eigenfunction. The resulting streamwise vorticity is shown in figure 14(b). Comparing the superposition of the wavepacket modes with the bi-global eigenfunction in figure 14(c), we see that the wavepackets capture the location and direction of the global mode, but lack the curvature of the wavefronts in the global structure near $y<0$ , $z\approx 0$ . A slice across the vortex core along $z=0$ is shown in figure 14(d), with the blue line representing the wavepacket results and the black line designating the global results. The red lines indicate where the individual wavepackets are centred (i.e. the value of $y_{\ast }$ ). A reasonable match between the bi-global analysis and the wavepacket analysis is observed.

Figure 14. Comparison of wavepackets to bi-global spectrum for the four wavepackets nearest to the eigenvalue $\unicode[STIX]{x1D714}=5.220-\text{i}0.501$ . (a) Shows the eigenvalue spectrum with the eigenvalue examined circled in red. (b) Shows the weighted sum over the four wavepackets from the wavepacket analysis. (c) Shows the bi-global eigenmode corresponding to the circled eigenvalue. (d) Is the projection of the disturbance velocity, $v$ , along the line $z=0$ , showing how the wavepackets (blue) compares with the bi-global spectrum (black). The red lines denote the wavepacket locations.

Due to the nature of wavepacket solutions and our access to an (approximate) analytical expression for the dispersion relation, we are in a position to explore and quantify the travel direction of modal wavepacket structure. We simply examine the group velocity, $\unicode[STIX]{x2202}\unicode[STIX]{x1D714}/\unicode[STIX]{x2202}\boldsymbol{k}$ , where $\boldsymbol{k}$ denotes the wavenumber vector with $\unicode[STIX]{x1D6FD}_{\ast }$ and $\unicode[STIX]{x1D6FE}_{\ast }$ aligning with the $y$ - and $z$ -directions, respectively. More specifically, an approximation of the group velocity can be extracted directly from the terms in the twist condition, $\unicode[STIX]{x2202}f/\unicode[STIX]{x2202}\unicode[STIX]{x1D6FD}_{\ast }$ and $\unicode[STIX]{x2202}f/\unicode[STIX]{x2202}\unicode[STIX]{x1D6FE}_{\ast }$ in (C 4). The group velocity of the dominant wavepacket (with the largest coefficient in the above superposition) shows that the disturbance travels inbound towards the vortex core, as indicated by arrows in figure 15(a,b). This inboard direction is congruent with the results from the bi-global analysis (see e.g. movie 4 of the supplemental material). Interestingly, although the velocity field appears outside the region of shear, the associated pressure field remains confined within the shear region as shown in figure 15(b). This may suggest a particular receptive behaviour whereby velocity disturbances in the far field trigger disturbances in the shear region via pressure perturbations, followed by a transfer of energy to the unstable discrete branch or transient growth. The interaction of free-stream disturbances with regions of shear is a common feature in receptivity studies (Saric et al. Reference Saric, Reed and Kerschen2002). The underlying mechanism behind this transfer from vorticity to pressure requires further research but may be a result of non-normality of the eigenvectors.

Figure 15. (a) The streamwise vorticity for the eigenvalue $\unicode[STIX]{x1D714}=5.220-\text{i}0.501$ and (b) the corresponding disturbance pressure. The solid contour lines denoting the location of the wake and the base of the arrow is the location of the wavepacket. The direction of the arrow indicates the group velocity, $\unicode[STIX]{x2202}\unicode[STIX]{x1D714}/\unicode[STIX]{x2202}\boldsymbol{k}$ , of the wavepacket.

Although figures 14 and 15 provide a single wavepacket, how the wavepackets vary across the continuous branch wedge is now examined qualitatively for the sake of brevity. When approaching the region near the analytical continuous branch in figure 13, the envelope of the wavepacket analysis broadens in the free stream, approaching the free-stream oscillation limit. Conversely, as the discrete branch is approached, the energy of the wavepacket is more compact and approaches the regions of shear. This physically aligns with continuous branch and discrete branch limits of the wavepacket analysis, matching well with the bi-global analysis and further providing confidence in the method.

This concludes the temporal stability analysis, which examined and identifies three distinct types of instabilities from the discrete spectrum, and assessed the continuous spectrum far from the shear flow. We then explored approximate wavepacket solutions. We now proceed to complement the temporal results with a spatial stability analysis.

4.1 Frequency parameterization

By sweeping the frequency from $\unicode[STIX]{x1D714}=2$ to 8 in increments of $\unicode[STIX]{x0394}\unicode[STIX]{x1D714}=0.5$ , we observe a number of instabilities including the addition of a vortex instability that contains high growth at low frequencies as shown in figure 16. Aside from this instability, the frequency sweep in figure 16 follows similar trends to the wavenumber sweep of the temporal analysis, shown in figure 4, containing principal wake, wake and vortex instabilities. As the spatial and temporal parametric sweeps contain similar structures, we choose to examine the frequency counterpart to the temporal analysis, which corresponds to $\unicode[STIX]{x1D714}=5.5$ indicated by the dashed line in figure 16.

Figure 16. Spatial growth rate, $-\unicode[STIX]{x1D6FC}_{i}$ , at $x=3$ with varying frequency $\unicode[STIX]{x1D714}$ , with $N_{z}=N_{y}=60$ . The four instabilities are labelled as wake mode, primary vortex mode, wake mode and secondary vortex mode according to their modal features. The dashed line indicates the wavenumber that we examined throughout the temporal analysis.

4.2 Eigenvalue spectrum and its comparison with the temporal spectrum

Due to the ansatz in (2.2), the eigenvalue spectrum for the spatial analysis inverts such that modally unstable eigenvalues reside in the lower half-plane $\unicode[STIX]{x1D6FC}_{i}<0$ , as shown by the grey region in figure 17(a). The discrete spectrum now sweeps over an area upward from $\unicode[STIX]{x1D6FC}\approx \unicode[STIX]{x1D714}+\text{i}\unicode[STIX]{x1D714}^{2}/Re$ , with discrete branches to its right and a single discrete branch to its left. The discrete branches right of the continuous spectrum contain a vortex and wake branch, analogous to their temporal counterparts. An additional third branch contains vortex modes to the farther right than the wake branch. The branch left of the continuous spectrum is analogous to the discrete branch to the right of the temporal continuous branch, which corresponds to higher-order azimuthal modes.

Figure 17. The spatial eigenvalue spectrum for $\unicode[STIX]{x1D714}=5.5$ in the complex $\unicode[STIX]{x1D6FC}$ -plane, (a), and in the complex phase-speed plane, (b). In (b), the equivalent temporal spectrum is shown for $\unicode[STIX]{x1D6FC}=5.5$ , denoted by the red circles. For comparison, the pertinent discrete branches are labelled.

The temporal and spatial spectra possess many similar features. For a more direct comparison, figure 17(b) compares the temporal spectrum (red circles) to the spatial spectrum (black dots) in the complex phase-speed plane, $c=\unicode[STIX]{x1D714}/\unicode[STIX]{x1D6FC}$ . In the $c$ -plane, the two spectra align, with the wake and vortex branches in the temporal analysis aligning with the corresponding branches in the spatial spectrum. The leftmost vortex branch in figure 17(b) aligns with stable modes in the temporal spectrum. The continuous branches sweep over approximately the same area, except that the spatial spectrum bends left as $c_{r}$ decreases – a phenomenon that is explained in § 4.4.

The alignment of the temporal and spatial eigenvalue spectra in $c$ -space has relevance to Gaster’s transformation (Gaster Reference Gaster1962), which holds for regions near the neutral stability curve. Gaster (Reference Gaster1962) showed that for small growth rates, the temporal and spatial growth rates can be related through the group velocity, $c_{g}=\unicode[STIX]{x2202}\unicode[STIX]{x1D714}_{r}/\unicode[STIX]{x2202}\unicode[STIX]{x1D6FC}_{r}$ . From the relatively good alignment of the temporal and spatial spectra in $c$ -space, we conclude that the deviation of the spectra is a result of the fairly large growth rates. Furthermore, the group velocity of the modes can be approximated, showing that these instabilities (and modes in general) are convective in nature.

4.3 Discrete spatial bi-global modes

The striking similarities between the temporal and spatial spectra suggest that the two spectra contain analogous eigenmodes on the matching branches. Thus, we expect the most unstable eigenvalue to represent a wake mode, similar to its temporal counterpart, which is confirmed in figure 18 (and movie 6 of the supplemental material). Comparing figure 18 with its temporal counterpart (figure 6), we observe that the disturbances remain localized to the wake region and contain nearly identical structures, which instead grow with downstream progression because of their instability. Furthermore, this instability has no spanwise zero crossings until the vortex is reached, again highlighting that the wake aft of the wing is the cause of this instability.

Figure 18. Bi-global spatial principal wake instability mode. Variables are plotted following the same format as in figure 6 with $Q=\pm 0.0001$ . For a movie of (c), see movie 6 in the supplemental material.

From these similarities, we expect that the modes farther down the wake branch contain higher-order spanwise structures, seen in the third instability shown in figure 19 (supplementary movie 7). This mode, lying close to its analogous temporal mode (see figure 7), contains a higher number of zero crossings, indicating a higher spanwise wavenumber. The same mode, however, decays more rapidly in the vicinity of the vortex. If we continue down this branch, we follow the same trend shown in the temporal section of higher-order spanwise structures in the wake region. Moreover, we argue that the reduced structure size and increased viscous effects cause the stabilization of these modes.

Figure 19. Bi-global spatial wake instability mode. Variables are plotted following the same format as in figure 6 with $Q=\pm 0.0001$ . For a movie of (c), see movie 7 in the supplemental material.

Continuing with the temporal analogues, along the vortex branch we observe that the disturbances are localized to the vortex region, with a wavepacket below it; compare figure 20 (supplementary movie 8) with figure 8. Contrary to the temporal mode, the vortex mode from the spatial analysis contains a complex vortex structure rather than the two-lobed structure that co-rotates with the base flow as in its temporal counterpart. The complex structure and wake aligns better with higher-order temporal vortex modes as shown in figure 9. Continuing down this branch confirms that the wavepacket contains higher-order modes.

Figure 20. Bi-global secondary spatial vortex instability mode. Variables are plotted following the same format as in figure 6 with $Q=\pm 0.0001$ . For a movie of (c), see movie 8 in the supplemental material.

From the spatial analysis, a new instability emerged, which is linked solely to the vortex region, with only fringe amounts of the disturbance localized in the wake, shown in figure 21 (supplementary movie 9). This instability contains a two-lobed structure where the inboard lobe originates from the region where the vortex and wake pinch off to separate, and the other lobe coming from the vortex core region. The rotation of this mode is difficult to discern with the disturbances traveling in the upward and inboard direction. Higher-order modes further down this branch contain smaller-scale structures in the transverse direction and comprise more disturbances in the wake region.

Figure 21. Bi-global primary spatial vortex instability mode. Variables are plotted following the same format as in figure 6 with $Q=\pm 0.0001$ . For a movie of (c), see movie 9 in the supplemental material.

The final mode we examine lies on a discrete branch to the right of the continuous spectrum in the $c$ -plane. This branch, termed the azimuthal branch in § 3.2, contains vortex modes with higher azimuthal wavenumbers. The mode chosen in figure 22 is an $(m=3)$ -mode that co-rotates with the flow (see supplementary movie 10) and is confined to the vortex region of the base flow. Modes on this branch may be significant at different frequencies, as Paredes (Reference Paredes2014) shows that higher-order modes contain higher growth rates in the spatial analysis of the Batchelor vortex; however, higher frequencies are beyond the scope of this study.

Figure 22. Bi-global spatial higher-order azimuthal mode. Variables are plotted following the same format as in figure 6 with $Q=\pm 0.0001$ . For a movie of (c), see movie 10 in the supplemental material.

4.4 Spatial continuous spectrum

The spatial continuous spectrum contains a peculiar bending when plotted in the complex $c$ -plane, shown in figure 17(b). Following an approach similar to Grosch & Salwen (Reference Grosch and Salwen1978), we examine the asymptotic limit of the spatial continuous spectrum (for small $\unicode[STIX]{x1D714}$ and $1/Re$ ), which in the complex $\unicode[STIX]{x1D6FC}$ -plane is parameterized by $\unicode[STIX]{x1D701}\in \mathbb{R}$ such that

(4.1) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6FC}=\unicode[STIX]{x1D714}+\text{i}\frac{1}{Re}(\unicode[STIX]{x1D714}^{2}+\unicode[STIX]{x1D701}). & & \displaystyle\end{eqnarray}$$

This is analogous to the temporal spectrum; however, differences become apparent after transforming the spectrum into the complex $c$ -plane. For the temporal case, we scale the spectrum by the inverse of the real wavenumber, $1/\unicode[STIX]{x1D6FC}$ . Conversely, in the spatial spectrum, equation (4.1) appears in the denominator. Finding the equivalent $c=\unicode[STIX]{x1D714}/\unicode[STIX]{x1D6FC}$ results in

(4.2) $$\begin{eqnarray}\displaystyle c=\frac{\unicode[STIX]{x1D714}^{2}-\text{i}\displaystyle \frac{\unicode[STIX]{x1D714}}{Re}(\unicode[STIX]{x1D714}^{2}+\unicode[STIX]{x1D701})}{\unicode[STIX]{x1D714}^{2}+\displaystyle \frac{(\unicode[STIX]{x1D714}^{2}+\unicode[STIX]{x1D701})^{2}}{Re^{2}}}, & & \displaystyle\end{eqnarray}$$

which corresponds to the black curve in figure 23 that represents the spatial continuous branch in the complex $c$ -plane. The curvature of this line matches the curvature of the spatial spectrum, illustrated by the dots in figure 23. Note that we solved for 5000 eigenvalues, which explains why the results from our analysis only partially cover the continuous branch. From these findings, we therefore conclude that the bending is a result of the inverse relation between $c$ and $\unicode[STIX]{x1D6FC}$ .

Figure 23. The eigenvalues from the bi-global analysis (red dots), the matched Batchelor vortex (blue circles), and the matched wake profile (yellow plus marks). The black line provides the analytical continuous branch given by (4.2) (black line).

4.5 Comparison of direct numerical simulation (DNS) stability results with canonical base flows

The stability analysis performed on the computationally obtained base flow imparts insight into the importance of the wake–vortex interaction on the stability of the intermediate field. We now provide a basis for comparison by juxtaposing the spatial spectrum from the computational base flow to matched canonical base flows; in particular the Batchelor vortex and wake velocity profile.

To match the base-flow field of the Batchelor vortex, we first azimuthally average the azimuthal and axial velocity components at $x=3$ and fit a Batchelor vortex to the average values by minimizing the least-squares residual. The result is the matched parameters of the Batchelor vortex, with the swirl strength $\unicode[STIX]{x1D705}=0.013$ , core radius $\unicode[STIX]{x1D6FF}=0.138$ , axial velocity deficit $\unicode[STIX]{x1D6FE}=0.229$ , which yields a swirl parameter $q=0.419$ . The axisymmetric base flow is then provided as

(4.3a,b ) $$\begin{eqnarray}\displaystyle V_{\unicode[STIX]{x1D703}}=\frac{\unicode[STIX]{x1D705}}{r}\left[1-\exp \left(-\frac{r^{2}}{\unicode[STIX]{x1D6FF}^{2}}\right)\right],\quad U=1-\unicode[STIX]{x1D6FE}\exp \left(-\frac{r^{2}}{\unicode[STIX]{x1D6FF}^{2}}\right), & & \displaystyle\end{eqnarray}$$

with the swirl parameter defined as $q=\unicode[STIX]{x1D705}/(\unicode[STIX]{x1D6FE}\unicode[STIX]{x1D6FF})$ . Figure 24 $(a,b)$ shows the comparison between the azimuthally averaged base flow (solid line) and the matched canonical fit (dashed line) for the azimuthal and axial velocity, respectively.

For the wake velocity profile, we extract the symmetry plane of the computed base flow and fit the profile given by Mattingly & Criminale (Reference Mattingly and Criminale1972) in a similar manner to the Batchelor vortex. The profile is

(4.4) $$\begin{eqnarray}\displaystyle U=U_{\infty }+(U_{c}-U_{\infty })\text{sech}^{2}(\unicode[STIX]{x1D70E}y), & & \displaystyle\end{eqnarray}$$

where $U_{\infty }=1.00$ , $U_{c}=0.626$ and $\unicode[STIX]{x1D70E}=7.84$ . The resulting fit is shown in figure 24(c). For more details of how the methodology of the fits, see Edstrand et al. (Reference Edstrand, Davis, Schmid, Taira and Cattafesta2016).

Figure 24. Comparison of base flow with fit canonical base flows. The solid line represents the DNS computed base-flow field and the dashed line denotes the canonical base flow. For panels (a,b), the DNS base flow was first azimuthally averaged prior to the Batchelor vortex being fit to the flow field; therefore, one must remember that the wake component averages out in the solid line plot in (a,b). In (c), we compare the symmetry plane of our computation to the base flow given by Mattingly & Criminale (Reference Mattingly and Criminale1972).

For the matched Batchelor vortex, we performed a bi-global spatial stability analysis on the base flow for the same parameters as the trailing vortex flow field. The resulting bi-global eigenvalue spectrum is shown in figure 23 denoted by the blue circles. Here, we ensured that the discrete branch was unaffected by the choice of the domain size. The spectrum shows two instabilities that aligns along the vortex branch and the azimuthal branch (see, e.g. figure 17 b) further reinforcing that the modes along these branches are vortex modes. When comparing the eigenfunctions in figure 25, we note that the disturbance of the computationally obtained base flow has strong wake–vortex interaction that is lacking in the Batchelor vortex instability. This lacking wake component, although expected as the base flow is axisymmetric, motivates the need for the wake and vortex component in the base flow; however a matched Batchelor vortex does result with a reasonable approximation of phase speeds and growth rates.

Figure 25. Comparison of the streamwise vorticity of the Batchelor vortex instability (a) with the bi-global vortex instability (b). Although the discrete branches align, the Batchelor vortex fails to capture the wake–vortex interaction.

For the wake velocity profile, the eigenvalue spectrum (denoted by the yellow plus marks in figure 23) aligns well with the wake branch of the computational base flow (again, see figure 17 b). The phase speed and growth rate of the principal wake instability, which contained uniform spanwise variation along the wake, aligns well with the canonical wake, indicating this is indeed a result of the wake aft of the wing. Furthermore, the eigenfunctions, not shown for brevity, match with the canonical value in both size and structure.

The comparisons of the computational base flow with the Batchelor vortex and the wake profile reinforces both the interpretation of the sub-branches in the discrete branch and the necessity to incorporate the wake and vortex structures in the base flow. We therefore conclude that although the canonical flow fields lack the necessary characteristics to fully illustrate the structure of the disturbance, the stability results do yield adequate approximations to phase speed and growth rates, indicating these canonical configurations are indeed a valuable first approximation.