2.1. Large eddy simulation

The database for both the data-driven and the theoretical mode calculations is obtained from a large eddy simulation performed with the commercial code STAR-CCM+ 14.02. A complete description of the numerical set-up can be found in Li et al. (Reference Li, Liang, Wang, Xiong, Chen, Yu and Chen2021b), here, only the essential information is presented.

A simplified version of the Intercity-Express 2, also known as the aerodynamic train model, is considered. The simulation is carried out on the scale of $1:10$, resulting in a height of the model of $H=0.36\,{\rm m}$, width of the model of $W=0.30\,{\rm m}$ and length of the model of $L=2.5\,{\rm m}$, as shown in figure 1(a). To simulate the relative motion between the train and the surrounding environment, the upstream boundary, located 10$H$ from the train head, is assigned the inflow velocity $U_\infty =4\,{\rm m}\,{\rm s}^{-1}$. Meanwhile, the ground boundary, with a distance to the train bottom surface of 0.15$H$, is defined with the same moving velocity. The resulting $\textit{Re}$ based on the height of the train and $U_\infty$ is $9.5\times 10^4$. The downstream boundary is located 30$H$ from the train tail, with a zero static pressure outlet condition. On the side and roof of the computational domain, the symmetry plane boundary condition is assigned, with a distance of 10$H$ from the train model. The coordinate system is shown in figure 1(a), with the origin located at the ground height of the train tail nose tip. The computational domain is then discretized using unstructured hexahedral volume meshes, with the wall-normal and wall-parallel distances expressed in viscous units, respectively, ${\Delta }y^+=0.16$ and ${\Delta }x^+={\Delta }z^+=28$ for cells attached to the train surface, as shown in figure 1(b). The total number of volume meshes used in the study is 46.8 million.

Figure 1. Large eddy simulation set-up. Reproduced from Li et al. (Reference Li, Liang, Wang, Xiong, Chen, Yu and Chen2021b). (a) Geometric model. (b) Distribution of computational grid. (c) Instantaneous snapshot of vortex structures; the arrow shows the flow direction.

In the current research, the large eddy simulation based on the wall-adapting local-eddy viscosity (WALE) subgrid-scale model is chosen. The use of a novel form of the velocity gradient tensor in the WALE subgrid-scale model allows for much more universal model coefficients compared with other widely used subgrid-scale models. Meanwhile, the WALE subgrid-scale model does not require any form of near-wall damping but automatically provides accurate scaling at the walls. More details about the WALE subgrid-scale model can be found in Nicoud & Ducros (Reference Nicoud and Ducros1999). An implicit unsteady segregated incompressible finite-volume solver is used, with the convective terms discretized based on a bounded central-differencing scheme, and the diffusion and turbulence terms discretized with the second-order upstream scheme. The time marching procedure is performed using the implicit second-order accurate three-time level scheme, with the discretized convective time step set to $0.0067t^*$ ($t^*=W/U_\infty$), which leads to a Courant–Friedrichs–Lewy number below unity in most of the computational grids. An instantaneous scene of the flow structures around the generic high-speed train, which briefly illustrates the formation of the turbulent wake, is shown in figure 1(c). Note that the LES results have been validated in Li et al. (Reference Li, Liang, Wang, Xiong, Chen, Yu and Chen2021b), where the simulated pressure coefficients are compared with experimental results. Here, to further enhance the confidence of the results presented in the paper, the LES data are further validated by a wind tunnel experiment, in terms of both the time-averaged flow velocity and turbulent statistics. The detailed comparisons are shown in Appendix A.

We started to collect the field data after the simulation had been run for $50t^*$. It can be seen from the time-history curves of global quantities shown in figure 2 that the flow has already reached the statistically stationary state at this moment. Then, the three-dimensional snapshot database was continuously collected from a square box extending from $x/W=-1.33$ upstream of the tail nose tip to $x/W=6.67$ downstream of the tail nose tip, $y/W=1.33$ from the centre plane of the train in both spanwise directions, and $z/W=1.33$ from the ground in the vertical direction, for the duration of $800t^*$. The time step between two consecutive snapshots is $0.1t^*$, resulting in a total of 8000 snapshots collected during the simulation. These values have been shown to be sufficient to produce well-converged SPOD results according to the parametric study shown in Li et al. (Reference Li, Chen, Liang, Liu and Xiong2021a).

Figure 2. Time histories of the aerodynamic force coefficients. (a) Drag force. (b) Lift force.

2.2. Mean flow properties

We consider the time-averaged field as the base state for the linear stability analysis and introduce it in this section. Furthermore, prior knowledge of the mean field will facilitate understanding of the physics associated with the extracted coherent structures. In figure 3(a), the time-averaged sectional streamlines in the wake are shown at several representative locations. Meanwhile, vortex regions in the wake are identified by the isosurface of $\varOmega$ (Liu et al. Reference Liu, Wang, Yang and Duan2016), coloured by the magnitude of the vorticity. The $\varOmega$-method works as a vortex identification criterion similar to the $Q$-criterion (Chong, Perry & Cantwell Reference Chong, Perry and Cantwell1990) and ${\lambda }_2$-criterion (Jeong & Hussain Reference Jeong and Hussain1995), however, it is less sensitive to threshold value to capture both strong and weak vortices in different cases. Note that the wake flow is symmetric about the central plane after long-term time averaging, so only the $y>0$ half is shown here for better visualization. In figure 3(b), to provide with more quantitative information, the time-averaged pressure coefficient and streamwise skin-friction coefficient along the central line of the train (also shown in this figure) are plotted. The $x$-axis is denoted by the distance from the starting point along the tail central line in the clockwise direction. The regions with negative streamwise skin-friction coefficients indicative of flow recirculation are highlighted with green patches.

Figure 3. Time-averaged mean flow distributions. (a) Flow structures around the train tail visualized by isosurface of vortex identification criterion $\varOmega =0.6$ and streamlines; the arrow shows the flow direction. (b) Distributions of pressure coefficient and streamwise skin-friction coefficient along the central line of the train tail in the clockwise direction; the greed patches highlight the regions with negative skin-friction coefficients.

The distribution of the time-averaged flow along the symmetry plane is illustrated first. It can be observed that, as the flow approaches the slanted surface of the tail, the adverse pressure gradient imposed by local flow acceleration forces the attached flow to separate from the surface at point $a$ (see figure 3b). However, the flow separated from point $a$ is highly deformed and fails to form a strong vortex region, since in figure 3(a), no obvious structure is identified at this location by the isosurface of $\varOmega$. Then, in figure 3(b), the re-attachment can be identified at point $b$, downstream of which the attached flow on the slanted surfaces gradually approaches the rear end. Further downstream, the strong adverse pressure near the tail nose point forces the attached flow on the slanted and bottom surfaces to respectively separate at point $c$ and point $d$, forming the spanwise vortex pair located just behind the tail nose point, as identified in figure 3(a).

In addition to the flow structures related to the separation across the symmetry plane, we can also observe in figure 3(a) that a pair of longitudinal vortices is located above the side edges of the tail, which is similar in nature to the C-pillar vortex in the Ahmed body wake flow (Zhang et al. Reference Zhang, Zhou and To2015; He et al. Reference He, Minelli, Wang, Dong, Gao and Krajnović2021a; Liu et al. Reference Liu, Zhang, Zhang and Zhou2021; Li et al. Reference Li, Cui, Jia, Li, Yang, Morzyński and Noack2022). This pair of longitudinal vortices is formed by flow separation from the side surface, which exerts a strong pressure-suction effect in this area and continuously rolls up flow from the slanted surface. As the longitudinal vortex propagates downstream, it gradually increases in diameter and lifts away from the tail surface toward the ground. Due to the strong downwash from the slanted surface, the trailing vortex is pushed away from the central symmetry plane as it travels downstream. From $x/W\approx 1.5$ onward, the longitudinal vortex structure attaches to the ground and then propagates nearly parallel to the ground in the downstream wake.

In general, the mean field is fully three-dimensional in the near-wake region, and its complexity gradually decreases downstream of the solid body, developing to be nearly parallel downstream of $x/W\approx 2.0$. The two main features, the spanwise recirculation bubble and the streamwise vortex pair, are likely to be related to the mean field instability due to the strong velocity gradient and will therefore be discussed further in the following content.

3.1. Spectral proper orthogonal decomposition

The SPOD approach is the frequency-domain counterpart of the standard POD approach. It seeks modes that optimally represent space–time flow statistics (Schmidt et al. Reference Schmidt, Towne, Rigas, Colonius and Brès2018; Towne et al. Reference Towne, Schmidt and Colonius2018). A brief overview of this method is provided in this part.

We denote the mean subtracted snapshots as

(3.1)\begin{equation} \boldsymbol{q}'_i=\boldsymbol{q}'(t_i)=\begin{bmatrix} u'(x,y,z,t_i) \\ v'(x,y,z,t_i) \\ w'(x,y,z,t_i) \\ p'(x,y,z,t_i)\end{bmatrix}, \end{equation}

where $i$ represents the number of snapshots, u, v, w and p are the normalized velocity components and pressure. To estimate the spectral contents using discrete Fourier transform (DFT), the snapshot database is first segmented into $n_{{blk}}$ overlapping blocks, with $n_{{DFT}}$ snapshots in each individual block. Here, $n_{{DFT}}=256$, along with an overlap of $50\,\%$, is used in our study, and the resulting angular frequency resolution is $\Delta \omega =0.245$ ($\omega$ is the angular frequency normalized by the factor of $U_{\infty }/W$). With the above parameters, $n_{{blk}}$ can be determined with the value of 61, which is sufficient for well-converged SPOD energy spectra and modes (Schmidt & Colonius Reference Schmidt and Colonius2020). The blocks are then Fourier transformed, and all Fourier realizations at the $k$th frequency are collected into a new data matrix

(3.2)\begin{equation} \boldsymbol{\hat{Q}}_k=\begin{bmatrix} \boldsymbol{\hat{q}}^{(1)}_k & \ \boldsymbol{\hat{q}}^{(2)}_k & \cdots & \boldsymbol{\hat{q}}^{(n_{{blk}})}_k \end{bmatrix}. \end{equation}

At this point, the orthogonal basis can be obtained by solving the eigenvalue problem defined by

(3.3) \begin{equation} \left(\boldsymbol{\hat{Q}}_{k}^H\boldsymbol{W}{\boldsymbol{\hat{Q}}_k}\right) / n_{{blk}}=\boldsymbol{U}_k\boldsymbol{\varLambda}_k\boldsymbol{U}_k^H, \end{equation}

where $\boldsymbol{\varLambda}_k$ and $\boldsymbol{U}_k$ are respectively the eigenvalue and eigenvector of this problem, and $\boldsymbol {W}$ is the diagonal matrix containing weight information of each flow quantity at each sampling node. Therefore, the weight matrix $\boldsymbol{W}$ defines the energy norm used in SPOD, thus determining the physical process to be highlighted (Colonius et al. Reference Colonius, Rowley, Freund and Murray2002). Here, the weight matrix is given as

(3.4)\begin{equation} \boldsymbol{W}=\int_{\varOmega}\begin{bmatrix} 1 & & & \\ & 1 & & \\ & & 1 & \\ & & & 0 \end{bmatrix}\,{{\rm d}}V \end{equation}

so that the turbulent kinetic energy (TKE) norm can be defined. The matrices $\boldsymbol {\varLambda }_k={\rm {diag}}(\lambda _{k_1},\lambda _{k_2},\ldots,\lambda _{k_{n_{blk}}},)$ contain the SPOD energies which are therefore based only on the TKE. Then both the velocity modes, as well as the associated pressure modes, can be recovered by

(3.5)\begin{equation} \boldsymbol{\varPhi}_k=\frac{1}{\sqrt{n_{{blk}}}}\boldsymbol{\hat{Q}}_k\boldsymbol{U}_k\boldsymbol{\varLambda}_k^{{-}1/2}. \end{equation}

In addition, since the described flow problem is subjected to the spanwise symmetry, fluctuations of the sampled flow field can be divided into symmetric and antisymmetric contributions (Hack & Schmidt Reference Hack and Schmidt2021). To properly isolate coherent structures defined by the two different types of fluctuations, an additive decomposition was applied to the collected snapshot data before extracting empirical modes. The symmetric contribution of the $i$th snapshot is defined as

(3.6)\begin{equation} \boldsymbol{q}'_S(t_i)=\frac{1}{2}\begin{bmatrix} u'(x,y,z,t_i)+u'(x,-y,z,t_i) \\ v'(x,y,z,t_i)-v'(x,-y,z,t_i) \\ w'(x,y,z,t_i)+w'(x,-y,z,t_i) \\ p'(x,y,z,t_i)+p'(x,-y,z,t_i) \end{bmatrix} \end{equation}

and the antisymmetric contributions is defined as

(3.7)\begin{equation} \boldsymbol{q}'_A(t_i)=\frac{1}{2}\begin{bmatrix} u'(x,y,z,t_i)-u'(x,-y,z,t_i) \\ v'(x,y,z,t_i)+v'(x,-y,z,t_i) \\ w'(x,y,z,t_i)-w'(x,-y,z,t_i) \\ p'(x,y,z,t_i)-p'(x,-y,z,t_i) \end{bmatrix} \end{equation}

A visualization of one snapshot of the fluctuating streamwise velocity field, together with the contributions from the symmetric and antisymmetric components, is shown in figure 4. The spanwise symmetrical and anti-symmetrical contributions of all collected samples are arranged into two separate data matrices. They are then independently solved following the procedures described above, to extract and analyse, respectively, the symmetric and antisymmetric empirical modes. Note that, due to the zero-integral property of an even and an odd function in the sampling domain, the symmetric and antisymmetric modes yield mutual orthogonality (Hack & Schmidt Reference Hack and Schmidt2021).

Figure 4. Symmetric and antisymmetric decomposition of a single snapshot. (a) Original field. (b) Symmetric component. (c) Antisymmetric component.

3.2. Spectral proper orthogonal decomposition energy spectra and modes

Spectral proper orthogonal decomposition solves the eigenvalue problems at each discrete frequency independently, and produces energy-ranked eigenvalues. Therefore, the energy distributions of different modes at each frequency can be best visualized in the form of a spectrum (Schmidt et al. Reference Schmidt, Towne, Rigas, Colonius and Brès2018). Figure 5(a,c) shows SPOD spectra for both symmetric and antisymmetric components. Meanwhile, the cumulative energy content and the percentage of energy accounted for by each mode as a function of frequency are shown in the right column. The symmetric component displays a higher overall energy level compared with the antisymmetric component, indicating its dominant role in the dynamics of the turbulent wake. The low-rank behaviour, characterized by a large separation between the first and second modes, also appears to be more pronounced in the symmetric component. In particular, in the angular frequency range of $2<\omega <4$ of the symmetric component, the first modes contribute more than $50\,\%$ of the fluctuating energy according to figure 5(b). The angular frequency of the dominant symmetric coherent structure is $\omega =3.437$, as shown in figure 5(a). The corresponding Strouhal number is $St=0.547$, based on the free-stream velocity $U_{\infty }$ and the train width $W$. Additionally, the spectral energy concentration at $\omega \approx 6.9$, and a less pronounced peak around $\omega \approx 1.7$, can be identified respectively. These two angular frequencies could be associated with the first harmonic and subharmonic of the fundamental mode. This will be further investigated by checking the mode bi-spectrum in § 3.3.

Figure 5. Spectral proper orthogonal decomposition mode energy spectrum. Panels (a,c) show the spectral curves with the shading of the line colour representing the increase in the mode number. Panels (b,d) show the percentage of energy that each mode accounts for as a function of frequency, with solid lines indicating that the cumulative energy represents $50\,\%$ and $90\,\%$ of the total energy at each frequency. (a,b) Symmetric component. (c,d) Antisymmetric component.

In figure 6, the spatial distributions of the leading symmetric SPOD modes at $\omega =3.437$, $\omega =1.718$ and $\omega =6.874$ are visualized based on the isosurfaces of the streamwise velocity components. A common feature of these modes is that they all behave very differently between the near- and far-wake regions, due to the complexity of the mean flow. For the leading SPOD mode at $\omega =3.437$, coherent vortex shedding is observed from the recirculation region just behind the train tail. The generated coherent structures move downstream and gradually approach the bottom boundary due to effect of the moving ground (Wang et al. Reference Wang, Minelli, Cafiero, Iuso, He, Basara, Gao and Krajnović2023). Once fully attached, they vanish and separate at the symmetry plane to evolve into far-wake coherent structures with nearly constant streamwise wavelength. For the leading mode at $\omega =1.718$, both alternate shedding of the spanwise vortex and co-shedding from the side surfaces are observed in the near wake. Similar to the most dominant SPOD mode, these structures slowly evolve into streamwise wavepackets when propagating to the far wake. Then, for the leading mode at $\omega =6.874$, no obvious structures can be identified in the near wake, and the far wake is mainly dominated by tilted wavepackets.

Figure 6. Spatial distributions of the leading symmetric SPOD modes visualized based on isosurface of the streamwise velocity components; (a) $\omega =3.437$, (b) $\omega =1.718$, (c) $\omega =6.874$.

Another important observation is that, with the increase of the frequency from $\omega =1.718$ to $\omega =6.874$, the spatial scales of the coherent structures gradually decrease. This phenomenon is more pronounced in the far wake, where the mean flow is nearly parallel and the coherent structures at different frequencies are all characterized by nearly constant streamwise wavenumbers. In particular, the streamwise wavelengths of the far-wake coherent structures at $\omega =1.718$ and $\omega =6.874$ can be observed to be approximately twice and half that of the most dominant SPOD mode, respectively. This phenomenon indicates a linear dispersion relation of the wake flow, which will be further discussed and confirmed in the following content.

To better quantify the frequency–wavenumber characteristics, a streamwise Fourier transform is applied to the streamwise velocity component of the SPOD modes, to convert the signal into the domain of streamwise wavenumber $\alpha$ (normalized by a factor of $1/W$). Due to the dominant role of the symmetric perturbation in the wake, only symmetric modes are considered. First, the power spectral density (PSD) of the streamwise velocity component of the leading mode is averaged along the $z$-axis following the expression

(3.8)\begin{equation} \bar{P}_1(\alpha,\omega,y)=\frac{\int_{z_{{min}}}^{z_{{max}}} \left|\hat{\varPhi}_u(\alpha,\omega,y,z)\right|\,{{\rm d}}z}{z_{{max}}-z_{{min}}}, \end{equation}

where $\hat {\varPhi }_u$ denotes the streamwise Fourier transform of the streamwise velocity mode, and $0\leq {z/W}\leq 1.3$. Meanwhile, we isolate the mechanism in near-wake and far-wake regions by using two different spatial window functions (Schmidt et al. Reference Schmidt, Towne, Rigas, Colonius and Brès2018) that constrain the signal to $0< x/W<1$ and $1< x/W<6.67$ when the Fourier transform is performed. The results of $\bar {P}_1$ at $\omega =3.437$ using the two different window functions are then shown in figures 7(a) and 7(b), respectively. In the near-wake region, the maximum PSD is located in the symmetric plane, with $\alpha \approx 8$. We can also observe that the high PSD value occurs over a wide range of $\alpha$, which can be attributed to the growth of the wavelength of the spanwise vortex shedding mode. Note that the non-zero values along $\alpha =0$ in figure 7(a) are due to spectral leakage in the discrete Fourier transform. In the far-wake region, the coherent structure has a nearly constant streamwise wavenumber, with the maximum PSD located away from the symmetry plane. This is consistent with what is shown in figure 6.

Figure 7. Streamwise wavenumbers of the SPOD modes: (a) $0< x/W<1$ at $\omega =3.437$; (b) $1< x/W<6.67$ at $\omega =3.437$; (c) Angular frequency–streamwise wavenumber diagram.

Then for the first SPOD modes at all discrete frequency points, a window function including both near wake and far wake is used to compute $\bar {P}_1$, followed by averaging along the $y$-axis, which can be expressed as

(3.9)\begin{equation} \bar{P}_2(\alpha,\omega)=\frac{\int_{y_{{min}}}^{y_{{max}}}\bar{P}_1(\alpha,\omega,y)\,{{\rm d}}y}{y_{{max}}-y_{{min}}}, \end{equation}

with $0\leq {y/W}\leq 1.2$, to construct the frequency–wavenumber diagram shown in figure 7(c). In general, a constant phase velocity of 0.83 can be observed for perturbation waves of all discrete frequencies considered in the research, which confirms the linear dispersion relation of the wake. This value lies between the convective velocity of the streamwise vortex and the outer free stream. Such a phase velocity is typical of the Kelvin–Helmholtz convective instability found in free-shear flow (Schmidt et al. Reference Schmidt, Towne, Rigas, Colonius and Brès2018). Since the frequency–wavenumber diagram includes both near-wake and far-wake instability waves, the phase velocity tends to decrease at frequencies with a pronounced low rank due to the near-wake mode.

3.3. Triadic interactions

The dynamic behaviour of a nonlinear flow system is significantly influenced by resonance phenomena driven by different mechanisms (Tang et al. Reference Tang, Cheng, Tong, Lu and Zhao2017). Triadic resonance, resulting from the quadratic nonlinearity of the Navier–Stokes equations, can play an important role in the energy transfer process in turbulent wake flows. Through this mechanism, coherent structures associated with the tonal components of vortex shedding can triadically interact with themselves as well as the background turbulence. As a consequence, the energy of the limit-cycle oscillation saturated from a fixed point supercritical bifurcation (Sipp & Lebedev Reference Sipp and Lebedev2007) is redistributed over different scales. Therefore, turbulent wake flows often exhibit a mixed tonal–broadband dynamics, comprising an underlying broadband spectrum and tonal components associated with vortex shedding. The broadband coherent dynamics with several identifiable peaks in the SPOD power spectrum shown in figure 5 indicates such behaviour, and is therefore further discussed in this section.

Bispectral mode decomposition (BMD) (Schmidt Reference Schmidt2020) detects triadic interactions by their characteristic phase coupling between two spectral components at frequencies $\omega _1$ and $\omega _2$, and a third frequency at $\omega _3$, obeying the resonance condition $\omega _1\pm \omega _2\pm \omega _3=0$. Bispectral mode decomposition extends classical bispectral analysis to multidimensional data, thereby simultaneously facilitating the detection of nonlinear triadic interactions from the complex mode bispectrum $\lambda _1(\omega _1,\omega _2)$, as well as the identification of the associated coherent flow structures as bispectral modes. For details on the derivation and implementation, the reader is referred to Schmidt (Reference Schmidt2020).

The mode bispectrum is then computed using the same spectral estimation parameters as in the SPOD of § 3.1. Since the flow is subject to the spanwise symmetry, triadic interactions can be further categorized into three different types: self-interaction of the symmetric component, self-interaction of the antisymmetric component and interaction between the symmetric and the antisymmetric components. Note that the third frequency component resulting from the self-interactions of both the symmetric and the antisymmetric components is symmetric, while the third frequency component resulting from the mutual interaction of symmetric and antisymmetric components is antisymmetric. Here, we confirm that the mode bispectra of the triadic interaction between the symmetric and antisymmetric components are overall lower in intensity and cannot be identified with distinguishable peaks, therefore will not be further discussed.

Hence, in figure 8, only the mode bispectra of the symmetric and antisymmetric self-interactions are presented. Consistent with the SPOD spectra, the mode bispectra also appear to be broadband. Here, the self-interaction of the symmetric component shown in figure 8(a) is identified with several distinct peaks. We refer to $\omega _f$ as the fundamental frequency detected in SPOD. Conceptually, the triplet $(\omega _f,0,\omega _f)$ on the mode bispectrum can be regarded as the coherent perturbation driven by the mean flow. However, due to the finite sampling frequency, the zero-frequency bin contains unresolved low-frequency components. In this case, the mode bispectrum on the $\omega _1$-axis should not be interpreted (Nekkanti et al. Reference Nekkanti, Maia, Jordan, Heidt, Colonius and Schmidt2022). A local maximum of the distribution can be found that corresponds to the sum interaction of the fundamental mode with itself, generating the first harmonic at twice the fundamental frequency via the triad $(\omega _f,\omega _f,\omega _{2f})$. At the same time, the doublets $(\omega _{2f},-\omega _f)$ and $(\omega _{({1}/{2})f},\omega _{({1}/{2})f})$, which include the first harmonic and subharmonic, respectively, interact with $\omega _f$. The triad $(\omega _f,-\omega _{({1}/{2})f},\omega _{({1}/{2})f})$, analogously, indicates an interaction with the subharmonic frequency. The mean field distortion, indicated by the peak values along line $\omega _1=-\omega _2$, can be identified over a wide frequency range. Meanwhile, in figure 8(b), a weak triad at $(\omega _{({1}/{2})f},\omega _{({1}/{2})f})$ can be identified, representing the sum self-interaction of the antisymmetric subharmonic contributing to the fundamental frequency in the symmetric component.

Figure 8. Mode bispectra of triadic interactions in both the sum and difference regions. (a) Self-interaction of the symmetric component. (b) Self-interaction of the antisymmetric component.

Selected dynamically relevant bispectral modes are shown in figure 9. Rows in this figure represent the constant third frequencies of $\omega _3=\omega _{({1}/{2})f}$, $\omega _3=\omega _{f}$ and $\omega _3=\omega _{2f}$. An important observation is that the modes along the diagonals are similar in shape, and resemble the SPOD modes at corresponding frequencies, as shown in figure 6. This confirms that the spatial structures involved in or generated by nonlinear triadic interactions with strong phase correlations are also the most energetic coherent structures (Nekkanti et al. Reference Nekkanti, Nidhan, Schmidt and Sarkar2023). The analysis confirms the expected nonlinear dynamics that is characterized by a fundamental vortex shedding frequency and different interactions of the fundamental mode at $(\omega _f,0,\omega _f)$ and its sub- and super-harmonics, including the mean flow distortion. Also apparent and consistent with the SPOD analysis, the train wake behaves stochastically, with the uncertainty of the fundamental frequency $\omega _f$ leading to broad spectral and bispectral peaks. Nonetheless, the continuous evolution and phase coupling of these structures results in the low-rank behaviour observed over a wide range of frequencies, as shown in figure 5. In addition to shedding light on the different nonlinear interactions that cumulatively lead to peaks in the SPOD energy spectra, the BMD analysis also reveals the self-interaction of the antisymmetric component at the subharmonic frequency in the mode bispectrum in figure 8(b) that couples the antisymmetric to the symmetric component.

Figure 9. Spatial distribution of the bispectral modes resulting from significant triads, visualized based on isosurface of the streamwise velocity component. (a) Self-interaction of the symmetric component. (b) Self-interaction of the antisymmetric component.

4.1. Linearized operator and treatment of the non-parallel flow

Coherent structures can be described by the triple decomposition (Reynolds & Hussain Reference Reynolds and Hussain1972), which leads to a further decomposition of the fluctuating component into coherent and stochastic parts as

(4.1)\begin{equation} \boldsymbol{q}'(\boldsymbol{x},t)=\boldsymbol{\tilde{q}}(\boldsymbol{x},t)+\boldsymbol{q}''(\boldsymbol{x},t), \end{equation}

where $\boldsymbol {\tilde {q}}(\boldsymbol {x},t)$ is the coherent fluctuation part and $\boldsymbol {q}''(\boldsymbol {x},t)$ is the stochastic fluctuation part. This decomposition is substituted into both the momentum and continuity equations, and both are time averaged and phase averaged. Then, by subtracting the time-averaged set of equations from the phase-averaged set of equations, the equations governing the evolution of coherent structures can be written (Reynolds & Hussain Reference Reynolds and Hussain1972)

(4.2)\begin{gather} \frac{\partial{\boldsymbol{\tilde{u}}}}{\partial{t}}+(\boldsymbol{\tilde{u}}\boldsymbol{\cdot}\boldsymbol{\nabla})\boldsymbol{\bar{u}}+(\boldsymbol{\bar{u}}\boldsymbol{\cdot}\boldsymbol{\nabla})\boldsymbol{\tilde{u}}={-}\boldsymbol{\nabla}{\tilde{p}}+\boldsymbol{\nabla}\boldsymbol{\cdot}\left(\frac{1}{\textit{Re}}(\boldsymbol{\nabla}+\boldsymbol{\nabla}^{\rm{T}})\boldsymbol{\tilde{u}}\right)-\boldsymbol{\nabla}\boldsymbol{\cdot}\left(\tau_R+\tau_N\right) \end{gather}

(4.3)\begin{gather} \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{\tilde{u}}=0. \end{gather}

Here, $\tau _N$ describes the quadratic interactions of the coherent perturbation. Considering that the energy contribution from this process is higher order, this term is neglected in the following. The term $\tau _R$ is the fluctuation of the stochastic Reynolds stresses related to the stochastic–coherent interaction, which contributes at leading order, according to the energy considerations of Reynolds & Hussain (Reference Reynolds and Hussain1972), and is therefore retained in the equation. However, this term is not known a priori and needs to be properly modelled to close the governing equation. In this paper, we use Boussinesq's eddy viscosity model as the closure. The Reynolds stresses are then expressed as

(4.4)\begin{equation} \tau_R={-}\nu_t\left(\boldsymbol{\nabla}+\boldsymbol{\nabla}^{\rm{T}}\right)\boldsymbol{\tilde{u}}. \end{equation}

Here, $\nu _t$ is the normalized eddy viscosity, which can be calculated using quantities of the LES mean flow. Since this approach yields an eddy viscosity for each independent Reynolds Stress component, a reasonable compromise is to take a value of $\nu _t$ that minimizes the mean squared error from the six constitutive relations using

(4.5)\begin{equation} \nu_t=\frac{\left\langle-\overline{\boldsymbol{u}'\boldsymbol{u}'}+2/3k\boldsymbol{I},\bar{\boldsymbol{S}}\right\rangle_{{F}}}{2\left\langle\bar{\boldsymbol{S}},\bar{\boldsymbol{S}}\right\rangle_{{F}}}, \end{equation}

with $\langle {\cdot },{\cdot }\rangle _{{F}}$ denoting the Frobenius inner product, $k$ the kinetic energy, $\boldsymbol {I}$ the identity matrix and $\bar {\boldsymbol {S}}$ the mean shear strain rate tensor. This approach has been widely used in the linear stability analysis of turbulent flows, as presented in Tammisola & Juniper (Reference Tammisola and Juniper2016), Rukes et al. (Reference Rukes, Sieber, Paschereit and Oberleithner2016), Kaiser et al. (Reference Kaiser, Poinsot and Oberleithner2018), Müller et al. (Reference Müller, Lückoff, Paredes, Theofilis and Oberleithner2020) and Kuhn, Soria & Oberleithner (Reference Kuhn, Soria and Oberleithner2021).

The linearized momentum and continuity equations for the coherent perturbation can be obtained as

(4.6)\begin{gather} \frac{\partial{\boldsymbol{\tilde{u}}}}{\partial{t}}+(\boldsymbol{\tilde{u}}\boldsymbol{\cdot}\boldsymbol{\nabla})\boldsymbol{\bar{u}}+(\boldsymbol{\bar{u}}\boldsymbol{\cdot}\boldsymbol{\nabla})\boldsymbol{\tilde{u}}={-}\boldsymbol{\nabla}{\tilde{p}}+\boldsymbol{\nabla}\boldsymbol{\cdot}\left(\left(\frac{1}{\textit{Re}}+\nu_t\right)(\boldsymbol{\nabla}+\boldsymbol{\nabla}^{\rm{T}})\boldsymbol{\tilde{u}}\right), \end{gather}

(4.7)\begin{gather} \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{\tilde{u}}=0. \end{gather}

These equations can be then rewritten as

(4.8)\begin{equation} \boldsymbol{\mathcal{L}}\boldsymbol{\tilde{q}}=0, \end{equation}

where $\boldsymbol {\mathcal {L}}$ is the operator of the linearized equations superimposing the spatial discretization and base state (Paredes et al. Reference Paredes, Hermanns, Le Clainche and Theofilis2013).

By assuming that the mean field has much smaller derivatives in the streamwise direction than in the transverse and vertical directions, the system can be Fourier transformed in the streamwise direction, following the streamwise weakly non-parallel flow assumption. The coherent perturbation can be then written as

(4.9)\begin{equation} \tilde{q}(x,y,z,t)=\hat{q}(y,z)\exp\left[{\rm{i}}\left(\alpha{x}-\omega{t}\right)\right]+{\rm c.c.}, \end{equation}

where $\boldsymbol {\hat {q}}$ is the complex eigenfunction, $\alpha =\alpha _{{r}}+\rm {i}\alpha _{{i}}$ is the complex streamwise wavenumber, $\omega =\omega _{{r}}+\rm {i}\omega _{{i}}$ is the complex angular frequency and c.c. is the complex conjugate. Substituting (4.9) into (4.8) enables stability analysis of the mean field in the cross-flow plane at different streamwise locations. The global stability characteristics can then be recovered based on the concept of absolute/convective instability and the global mode wavemaker (Huerre & Monkewitz Reference Huerre and Monkewitz1990).

However, using the local approach to predict the global mode has been reported to be less accurate than the direct global approach when the base flow is strongly non-parallel (Juniper et al. Reference Juniper, Tammisola and Lundell2011; Juniper & Pier Reference Juniper and Pier2015). In the current case, where a fully developed three-dimensional base flow is considered, the parallel flow assumption is likely to introduce uncertainty into the results. Therefore, in this work, we also make further treatment to approximate the non-parallelism of the three-dimensional base flow.

The parallel flow assumption is checked by visualizing the streamwise component of the leading SPOD mode on a horizontal plane, as shown in figure 10. Here, the coherent perturbations can be observed with clear wavepacket structures; however, they do not travel strictly in the streamwise direction but follow oblique paths in both the near and far wake. This is caused by the downwash flow from the slanted tail surface which gradually separates the wake structures as it propagates downstream, as can be seen in the mean flow structures shown in figure 3.

Figure 10. Two-dimensional slice of the streamwise component of three-dimensional SPOD mode showing the oblique downstream travelling wavepackets. The streamlines are drawn based on the vector field ($\cos \theta, \sin \theta$) at each computational node.

We intend to account for the obliqueness of the coherent structure in the linear modelling. To do this, the $x$-dependence of the eigenfunction $\boldsymbol {\hat {q}}$ has to be re-considered. Then, for a given location $(x_0,y_0,z_0)$, and a streamwise distance $\Delta {x}$ small enough, there exists a set of $(\Delta {y},\Delta {z})$ so that the eigenfunction fulfils

(4.10)\begin{equation} \hat{q}(x_0,y_0,z_0)=\hat{q}(x_0+\Delta{x},y_0+\Delta{y},z_0+\Delta{z}). \end{equation}

This set of $(\Delta {y},\Delta {z})$ is related to the obliqueness of the travelling coherent structure, and (4.10) can be replaced by

(4.11)\begin{equation} \hat{q}(x_0,y_0,z_0)=\hat{q}(x_0+\Delta{x},y_0+\tan\theta\Delta{x},z_0+\tan\gamma\Delta{x}), \end{equation}

where $\tan \theta$ and $\tan \gamma$ respectively represent the convection direction relative to the symmetric plane and the horizontal plane. By applying a first-order Taylor expansion to the right-hand side of (4.11), the left-hand side can be cancelled and the expression can be arranged into

(4.12)\begin{equation} \frac{\partial{\boldsymbol{\hat{q}}}}{\partial{x}}+\tan\boldsymbol{\theta}\frac{\partial{\boldsymbol{\hat{q}}}}{\partial{y}}+\tan\boldsymbol{\gamma}\frac{\partial{\boldsymbol{\hat{q}}}}{\partial{z}}=0. \end{equation}

Therefore, an $x$-derivative applied to the state vector $\boldsymbol {\hat {q}}$ is used to account for the oblique travelling direction.

To determine $\tan \theta$ and $\tan \gamma$, we replace $\boldsymbol {\hat {q}}$ with mean field quantities $\boldsymbol {\bar {q}}$ in (4.12), by assuming that the perturbation waves follow the mean flow convection. Note that, for each node in the computational domain, the four flow variables are used to construct the linear equation system for $\tan \theta$ and $\tan \gamma$ and the final results are obtained using the least-squares solution of the overdetermined linear equation system.

Finally, the convection direction calculated from the mean field is validated by comparison with the oblique path of the SPOD mode. In figure 10, streamlines based on the calculated local angle $\theta$ are visualized, by decomposing $\tan \theta$ at each computational node into streamwise ($\cos \theta$) and transverse ($\sin \theta$) components using trigonometric functions. It can be observed that the vector field agrees well with the travelling direction of the SPOD wavepackets. Therefore, this $x$-derivative is assumed to be reasonable to account for the obliqueness of the coherent structure in the linear modelling. Note that the results of the stability analysis without the non-parallelism modelling are also computed and presented in Appendix B for comparison.

4.2. Spatio-temporal stability approach

For the linear stability analysis, the linear operator $\boldsymbol {\mathcal {L}}$ is rearranged to construct the generalized eigenvalue problem. In this work, both the temporal and spatial stability formulation is needed. Therefore, either the temporal stability form

(4.13)\begin{equation} \boldsymbol{\mathcal{A}}(x,\alpha)\boldsymbol{\hat{q}}=\omega\boldsymbol{\mathcal{B}}(x)\boldsymbol{\hat{q}}, \end{equation}

or the spatial stability form

(4.14)\begin{equation} \boldsymbol{\mathcal{A}}(x,\omega)\boldsymbol{\hat{q}}=\alpha\boldsymbol{\mathcal{B}}(x)\boldsymbol{\hat{q}}, \end{equation}

is constructed.

In the temporal stability form (4.13), the streamwise wavenumber $\alpha$ is fixed to a real value, and the eigenvalue problem is solved for a complex $\omega$, the real part of which is the angular frequency and the imaginary part is the temporal growth/decay rate. On the contrary, in the spatial stability form (4.14), a real angular frequency $\omega$ is imposed to search for the complex $\alpha$, where the real part corresponds to the streamwise wavenumber, while the imaginary part is the spatial amplification/damping rate. Note that, in the spatial stability equation, quadratic terms appear with respect to $\alpha$. This problem is solved using the companion matrix method (Bridges & Morris Reference Bridges and Morris1984), which increases the size of the eigenvalue problem compared with the temporal analysis. The complete expression of the operators $\boldsymbol {\mathcal {A}}$ and $\boldsymbol {\mathcal {B}}$ for both temporal and spatial analysis can be found in Appendix C.

To determine the global stability of the mean flow from a local analysis, a so-called spatio-temporal analysis using (4.13) must be performed to distinguish between convective and absolute instability (Huerre & Monkewitz Reference Huerre and Monkewitz1990). In this case, both $\alpha$ and $\omega$ are complex valued. Conceptually, the flow is said to be stable if all perturbations decay in time throughout the entire domain after a localized impulse. Convectively unstable flow gives rise to perturbations that grow in time but convect away from the impulse location, so that the perturbations eventually decay to zero at each spatial location in the long-time limit. For an absolutely unstable flow, the perturbations grow both upstream and downstream of the impulse location, contaminating the whole spatial domain in the long-time limit.

Based on the previous definitions, the convective/absolute nature of a local velocity profile can be determined from the time-asymptotic behaviour of the perturbations that remain at the impulse location, that is, for perturbations with zero group velocity: $\partial \omega /\partial \alpha =0$, which is the definition of a saddle point in the complex $\alpha -$plane. Therefore, valid saddle points on the complex $\alpha$-plane that satisfy the Briggs–Bers pinch-point criterion (Briggs Reference Briggs1964) must be identified. Therefore, the flow is locally absolutely unstable if the absolute growth rate, given by the imaginary part of $\omega$ at the most unstable valid saddle point, is positive. If the absolute growth rate is negative, the flow is convectively unstable or stable (Huerre & Monkewitz Reference Huerre and Monkewitz1990; Rees & Juniper Reference Rees and Juniper2010; Juniper et al. Reference Juniper, Hanifi and Theofilis2014; Rukes et al. Reference Rukes, Sieber, Paschereit and Oberleithner2016; Kaiser et al. Reference Kaiser, Poinsot and Oberleithner2018; Demange, Chazot & Pinna Reference Demange, Chazot and Pinna2020).

In a spatially developing flow, a region of absolute instability is a necessary (but not sufficient) condition for global instability (Monkewitz et al. Reference Monkewitz, Huerre and Chomaz1993; Chomaz Reference Chomaz2005). To further link the local absolute instability to global instability, the absolute growth rate needs to be tracked along the streamwise direction to determine the wavemaker, the location where the global mode is selected. This method has been extensively used for one-dimensional local stability analysis in comparison with bi-global stability analysis (Giannetti & Luchini Reference Giannetti and Luchini2007; Juniper et al. Reference Juniper, Tammisola and Lundell2011; Juniper & Pier Reference Juniper and Pier2015; Kaiser et al. Reference Kaiser, Poinsot and Oberleithner2018).

4.3. Solving the eigenvalue problem

To solve the eigenvalue problem numerically, cross-flow planes with the dimensions of $0\leq y/H\leq 1.3$ and $0\leq z/H\leq 1.3$ are discretized using Chebyshev spectral collocation methods. This approach has been successively applied to linear stability analysis by Khorrami (Reference Khorrami1991), Parras & Fernandez-Feria (Reference Parras and Fernandez-Feria2007), Oberleithner et al. (Reference Oberleithner, Sieber, Nayeri, Paschereit, Petz, Hege, Noack and Wygnanski2011) and Demange et al. (Reference Demange, Chazot and Pinna2020). Detailed descriptions or practical guides to spectral collocation methods can be found in Khorrami, Malik & Ash (Reference Khorrami, Malik and Ash1989) and Trefethen (Reference Trefethen2000).

To reduce the numerical cost, we further exploit the symmetry of the mean field with respect to the vertical $x$–$z$ plane. This allows us to use only half of the wake plane instead of the full plane, to compute only transversely symmetric or antisymmetric eigenmodes when appropriate boundary conditions are applied (Zampogna & Boujo Reference Zampogna and Boujo2023). As shown in § 3, the symmetric perturbations are dominant and potentially related to a global instability, so only symmetric eigenmodes are considered. The corresponding boundary conditions are given as

(4.15a)$$\begin{gather} \frac{\partial{\boldsymbol{\hat{u}}}}{\partial{y}}=\frac{\partial{\boldsymbol{\hat{w}}}}{\partial{y}}=\frac{\partial{\boldsymbol{\hat{p}}}}{\partial{y}}=0,\boldsymbol{\hat{v}}=0, \quad {\rm on} \ y/W=0, \end{gather}$$

(4.15b)$$\begin{gather}\frac{\partial{\boldsymbol{\hat{u}}}}{\partial{y}}=\frac{\partial{\boldsymbol{\hat{v}}}}{\partial{y}}=\frac{\partial{\boldsymbol{\hat{w}}}}{\partial{y}}=\frac{\partial{\boldsymbol{\hat{p}}}}{\partial{y}}=0, \quad {\rm on}\ y/W=1.3, \end{gather}$$

(4.15c)$$\begin{gather}\boldsymbol{\hat{u}}=\boldsymbol{\hat{v}}=\boldsymbol{\hat{w}}=0,\frac{\partial{\boldsymbol{\hat{p}}}}{\partial{z}}=0, \quad {\rm on} \ z/W=0, \end{gather}$$

(4.15d)$$\begin{gather}\boldsymbol{\hat{u}}=\boldsymbol{\hat{v}}=\boldsymbol{\hat{w}}=\boldsymbol{\hat{p}}=0, \quad {\rm on} \ z/W=1.3. \end{gather}$$

Here, (4.15a) determines eigenmodes to be transversely symmetric. At the wall, homogeneous Dirichlet boundary conditions are imposed for the velocity components. For pressure, by substituting the homogeneous Dirichlet conditions for velocity into (4.6), a compatibility condition can be obtained (Theofilis, Duck & Owen Reference Theofilis, Duck and Owen2004). Here, assuming $\partial ^2{\boldsymbol {\hat {w}}}/\partial {z}^{2}=0$ gives a homogeneous Neumann condition for the pressure. For far-field boundaries, the upper boundary is set to a homogeneous Dirichlet condition, while the side boundary is set to a homogeneous Neumann condition which is necessary to predict far-wake eigenmodes.

The Krylov–Schur algorithm (Stewart Reference Stewart2002), which serves as an improvement on traditional Krylov subspace methods such as the Arnoldi and Lanczos algorithms, is used to obtain a subset of solutions to the eigenvalue problem. This requires an initial guess of the physical eigenvalue, which can be derived from the dispersion relation shown in § 3.2. To discard spurious eigenmodes caused by the discretization, two criteria are applied: first, all eigenmodes that do not diminish when approaching the upper boundary are discarded. Second, since spurious eigenvalues are very sensitive to discretization, a convergence study of eigenvalues computed using different grid resolutions is used as a criterion for filtering spurious eigenvalues.

The results of the convergence study based on temporal stability analysis are shown in figure 11. Two representative cross-flow planes, located in the near wake and far wake, respectively, are considered. Real streamwise wavenumbers of 10 and 5 are, respectively, imposed in the two eigenvalue problems, which are then discretized using two different grid resolutions and solved for the subset of the eigenvalue spectrum. Note that the real streamwise wavenumbers are chosen according to the wavelength distribution of the dominant SPOD mode shown in figure 7, so that the SPOD peak frequency can be set as the initial value of the Krylov–Schur iteration procedure. As shown in figure 11, both spectra feature continuous branches and a set of discrete modes. In general, continuous branches are made up of spurious eigenmodes caused by discretization and physical eigenmodes that are highly stable.

Figure 11. The spectrum of eigenvalues obtained with discretization points of 8100 ($\square$) and 10 000 ($\circ$, blue), the physical eigenvalues are marked with $\ast$, red. (a) Results based on cross-flow plane at $X/W=0.1$ and imposing $\alpha =10$. (b) Results based on cross-flow plane at $X/W=3.5$ and imposing $\alpha =5$.

It can be observed that the locations of these eigenmodes in continuous branches vary significantly when the resolution of the grid is changed. On the contrary, the locations of the discrete modes remain almost stationary with different discretization settings; hence, the discrete modes are considered as physical eigenmodes that can potentially contribute to global instability. In particular, the near-wake cross-flow plane features one unstable mode and two stable modes, while the far-wake plane features two stable modes. Note that these physical eigenvalues do not necessarily correspond between the near-wake and far-wake cross-flow planes. When tracking along the streamwise direction, the physical eigenvalues may become highly stable and fall into the spurious region, and new physical eigenvalues may emerge due to the complexity of the base flow.

4.4. Absolute/convective stability analysis

4.4.1. Streamwise evolution of temporal growth rate

Before considering the absolute/convective nature of the instability, it is necessary to know which part of the flow is temporally unstable. To this end, the complex frequencies of the previously identified physical modes at $x/H=0.1$ and $x/H=3.5$ are tracked when varying the real streamwise wavenumber. The exact maximum temporal growth rate of each mode is then determined based on a Newton–Raphson iterative method (Ypma Reference Ypma1995). Figure 12 shows the branches of the most unstable modes at the two streamwise locations, with the associated maxima marked by asterisks. These maxima are then tracked along the streamwise direction, with a spacing of $\Delta {x}/W=0.002$, by repeating the iteration procedure. As mentioned above, due to the complexity of the base flow, one physical eigenmode may occur only in a certain range of streamwise location. Therefore, several more cross-flow planes are used to compute the physical eigenvalues, and then the tracking process is repeated to account for the maximum temporal growth rates of all physical eigenmodes in the entire wake.

Figure 12. Frequency (a) and temporal growth rate (b) of the of the temporally most unstable modes at $x/H=0.1$ (dashed line) and $x/H=3.5$ (solid line) as a function of real wavenumber. The maximum temporal growth rates of the two eigenmodes are marked by asterisks.

In figure 13, the maximum temporal growth rates of these physical eigenmodes are displayed as a function of streamwise location. Only the unstable regime ($\omega _{{i}}>0$) is shown here. Three different unstable eigenvalue branches are identified in the entire wake, each dominating in different streamwise regions. Accordingly, the flow becomes temporally unstable very close to the tail of the train and remains unstable in the wake. Based on the spatial locations where these branches become unstable, we conceptually divide the wake into the near-, middle- and far-wake regimes and name these branches accordingly. It can be found in figure 13(b) that the near-wake eigenmode is temporally more unstable than all the downstream eigenmodes. According to the mean flow distribution shown in figure 13(a), this eigenmode branch is attributed to the transverse recirculation zone located right behind the tail of the train. The far-wake eigenmode can be observed across a wide streamwise distance and is located close to the extension of the streamwise vortex pair. This branch has the largest temporal growth rate at $x/W\sim 0.4$, and gradually stabilizes as it extends into the far wake; however, it remains still unstable. The middle-wake mode becomes unstable at $0.5< x/W<1.4$. At this location, with respect to the mean field, the flow from the slanted tail surface can be observed attached to the ground according to figure 13(a).

Figure 13. Maximum temporal growth rate as functions of streamwise location. (a) Mean flow distribution in the central plane; (b) maximum temporal growth rate of near-wake (------, red), middle-wake (------, green) and far-wake (------, blue) eigenmodes with real $\alpha$ imposed.

4.4.2. Spatio-temporal stability analysis in the near wake

Spatio-temporal analysis is performed to compute the contour of $\omega$ in the complex $\alpha$-plane, so as to find valid saddle points following the Briggs–Bers criterion. Since the branch in the near-wake region is temporally more unstable than the other branches (figure 13), we start with the near-wake branch.

Figure 14 shows the contour of the complex frequency in the complex $\alpha$-plane, at the streamwise location of $x/W=0.08$. At this location, the near-wake eigenmode reaches its maximum temporal growth rate. In this figure two saddle points are found, $S_0$ and $S_1$. However, valid saddles must be pinched between an $\alpha ^+$ and an $\alpha ^-$ branch (Huerre & Monkewitz Reference Huerre and Monkewitz1990). A simple rule based on the Briggs–Bers pinch-point criterion (Briggs Reference Briggs1964) can be applied by checking the isocontours of the spatio-temporal growth rate (figure 14b): starting from the saddle point, the contours of the growing $\omega _{{i}}$ must reach, respectively, the $\alpha _{{i}}>0$ and $\alpha _{{i}}<0$ half-planes. Here, the validity of both $S_0$ and $S_1$ can be confirmed, with $S_0$ representing the shorter travelling wave at higher frequency, while $S_1$ represents the longer travelling wave at lower frequency. In the long-time limit, the nature of the instability is determined by $S_0$, which has a significantly higher absolute growth rate than $S_1$. Then the absolute frequency at this location can be determined as $\omega _0=3.2904+0.0685{\rm {i}}$, which shows that the flow is absolutely unstable at this position.

Figure 14. Contour of the complex angular frequency in the complex $\alpha$-plane at $x/W=0.08$. Saddle points on the complex $\alpha$-plane are marked by a red circle. (a) The real component; (b) the imaginary component; the $\omega _{{i}}=0$ isocontour is highlighted in red.

Furthermore, to ensure that all regions with absolute instability have been taken into account, the absolute growth rate has been determined for all unstable branches at several streamwise locations. The near-wake branch was the only one to reveal absolute instability.

4.5. Determining the global mode wavemaker in the near wake

To determine the global mode, the saddle points are tracked in the streamwise direction to find the global wavemaker. This is done in the local analysis framework based on the frequency selection criterion. Several criteria have been evaluated in Pier (Reference Pier2002), showing that the criterion for a linear global mode introduced by Chomaz, Huerre & Redekopp (Reference Chomaz, Huerre and Redekopp1991) agrees best with nonlinear direct numerical simulations. This criterion is obtained from an analytical continuation of the absolute frequency curve into the complex $x$-plane. The wavemaker region is then represented by the saddle point on the complex $x$-plane (Huerre & Monkewitz Reference Huerre and Monkewitz1990) defined as

(4.16)\begin{equation} \frac{\partial\omega_0}{\partial{x}}(x_{{s}})=0. \end{equation}

The global complex angular frequency is then given by the value of the absolute frequency at this saddle point

(4.17)\begin{equation} \omega_{{g}}=\omega_0(x_{{s}}). \end{equation}

For this purpose, the saddle point $S_0$ found in figure 13 is tracked along the streamwise direction. The three-point Taylor series expansion algorithm (Rees Reference Rees2010) is used to approach $\partial \omega /\partial \alpha =0$ during the tracking process. Then the absolute frequency as a function of streamwise location is analytically continued in the complex $x$-plane using the Padé polynomial, which has been shown to be well behaved in the complex plane (Cooper & Crighton Reference Cooper and Crighton2000). The Padé polynomial takes the form of

(4.18)\begin{equation} f(x)=\frac{P(x)}{Q(x)}=\frac{a_0+a_1x+\cdots+a_nx^n}{1+b_1x+\cdots+b_mx^m}. \end{equation}

To determine the order of the polynomials used in the current work, the procedure described in Juniper et al. (Reference Juniper, Tammisola and Lundell2011) is followed. Polynomials of order 10 have been proven to be sufficient to give a converged approximation of the saddle point.

Figure 15 shows the main results of the spatio-temporal stability analysis. For reference, the mean field in the central plane is shown in figure 15(a). The following rows include the results of the procedures described above to determine the global instability using local analysis. In figure 15(b,c), the imaginary and real parts of the absolute complex angular frequency are shown as a function of streamwise location, respectively. A small region of absolute instability is detected, which is located in the recirculation region. The tenth-order Padé polynomials show acceptable agreement with the absolute frequency curve of the linear stability analysis. The extrapolated complex $x$-plane is then visualized in figure 15(d). The saddle point appears to be very close to the real axis, with $x_{{s}}=0.0733+0.0010{\rm {i}}$, indicating a wavemaker located within the region of absolute instability. The complex global angular frequency associated with this saddle point is $\omega _{{g}}=3.2702+0.0797{\rm {i}}$.

Figure 15. (a) Streamlines and through-plane vorticity of the mean field in the near-wake region. (b) Streamwise evolution of the absolute growth rate. (c) Streamwise evolution of the absolute angular frequency. (d) Contours of $\omega _{{0,i}}$ extruded in the complex $x$-plane, formed by the fitted Padé polynomial. The saddle point on the complex $x$-plane is marked by a red circle. (e) The imaginary component of local complex streamwise wavenumber of $\alpha ^+$ and $\alpha ^-$ branches calculated at the global frequency $\omega _g$. ( f) The same as (e) with the real component shown. The switching locations between $\alpha ^+$ and $\alpha ^-$ branches are marked by a red vertical line in each panel.

The theoretically predicted global mode frequency can be compared with the frequency of the most energetic SPOD mode, as presented in table 1. Meanwhile, the global frequency and wavemaker location predicted by the stability analysis without the non-parallelism approximation (see Appendix B) are also shown for comparison. In general, the global frequency predicted by the stability analysis without the non-parallelism approximation shows a larger deviation from the SPOD result, compared with that with the non-parallelism approximation. By including the spatial variation of the mean field into the stability analysis, the prediction on the global stability has been significantly improved. With the empirical value of $\omega _{{SPOD}}=3.437$, we find that the theoretical prediction deviates by an error of less than $5\,\%$. This is in very good agreement, given the uncertainties introduced by the non-parallelism and eddy viscosity modelling. Moreover, we expect the global mode to be marginally unstable, representing a limit-cycle oscillation when stability analysis is performed on a temporally averaged mean flow (Noack et al. Reference Noack, Afanasiev, Morzyński, Tadmor and Thiele2003; Barkley Reference Barkley2006; Oberleithner et al. Reference Oberleithner, Sieber, Nayeri, Paschereit, Petz, Hege, Noack and Wygnanski2011; Rukes et al. Reference Rukes, Sieber, Paschereit and Oberleithner2016; Tammisola & Juniper Reference Tammisola and Juniper2016; Kaiser et al. Reference Kaiser, Poinsot and Oberleithner2018). However, as reported in Giannetti & Luchini (Reference Giannetti and Luchini2007), Juniper et al. (Reference Juniper, Tammisola and Lundell2011) and Juniper & Pier (Reference Juniper and Pier2015), it is a common feature of local analysis to overpredict the growth rate of the linear global mode. In our case, with a growth rate value of 0.0797, the relative error with respect to the real frequency (3.2702) is less than $3\,\%$, therefore, it is reasonable to say that the mode is marginally unstable. Overall, considering the fact that the flow is not truly parallel and the intrinsic assumption of linearized mean field analysis, it can be concluded that the theoretical mode identifies the dominant SPOD mode as a global mode at limit cycle with remarkable accuracy.

Table 1. Global frequency and wavemaker location predicted by different approaches.

The global mode is formed by the switching between the $\alpha ^+$ and $\alpha ^-$ branches at the global frequency. Therefore, to further truly localize the global mode wavemaker, a spatial stability analysis (4.14) imposing $\omega =\omega _{{g}}$ is conducted to compute the $\alpha ^+$ and $\alpha ^-$ branches (Juniper & Pier Reference Juniper and Pier2015). Since the saddle point $x_{{s}}$ on the complex $x$-plane is very close to the real axis, the location of maximum structural sensitivity should also be close to $x_{{s,r}}$. Therefore, in practice, it can be quite straightforward to find the branch pairs by computing the eigenvalue problem at $x/W=x_{{s,r}}$, and following $\alpha _{{i}}^+-\alpha _{{i}}^-\approx 0$ in the streamwise direction.

In figure 15(e,f), the imaginary and real components of the $\alpha ^+$ and $\alpha ^-$ branches are presented, respectively. The switch between the two branches can be identified to take place at $x/W=0.0731$. The location of the global mode wavemaker is marked by a red vertical line in all panels of figure 15. This location is also within the recirculation region and is nearly identical to the location of the saddle point $x_{{s}}$. The direct global mode then follows the $\alpha ^-$ branch upstream of the wavemaker, and the $\alpha ^+$ branch downstream. On the contrary, the adjoint global mode follows the $\alpha ^+$ branch upstream of the wavemaker, and $\alpha ^-$ branch downstream (Juniper & Pier Reference Juniper and Pier2015). This result also indicates a spatially growing global mode within the recirculation region, with the growth rate gradually decaying to zero as it approaches the downstream boundary of this region.

In summary, we identify a global mode with a frequency very close to the dominant SPOD mode and a growth rate close to zero. This suggests that the dominant SPOD mode is a manifestation of a global mode at the limit cycle. The wavemaker of this mode is located in the recirculation zone very close to the tail tip. It acts as the origin for the entire coherent wavepacket that propagates far downstream, in a region where the flow is convectively unstable. The spatial shape and mechanisms of the global mode and comparison with the SPOD modes will be shown and discussed in § 5.

4.6. Structural sensitivity based on the adjoint method

To further analyse where and how intrinsic feedback causes global instability (Chomaz Reference Chomaz2005; Giannetti & Luchini Reference Giannetti and Luchini2007), we calculate the structural sensitivity using adjoint linear stability analysis. The structural sensitivity describes the sensitivity of the direct global mode to changes in the linearized Navier–Stokes operator, e.g. base flow modification (Marquet et al. Reference Marquet, Lombardi, Chomaz, Sipp and Jacquin2009). Therefore, this concept is critical for the development of efficient flow control strategies (Giannetti & Luchini Reference Giannetti and Luchini2007; Müller et al. Reference Müller, Lückoff, Paredes, Theofilis and Oberleithner2020).

In the global framework, structural sensitivity is proportional to the overlap between the direct and adjoint global modes (Chomaz Reference Chomaz2005). Following Giannetti & Luchini (Reference Giannetti and Luchini2007), the $L^2$ norm of the sensitivity tensor is equivalent to the expression

(4.19)\begin{equation} \lambda(\boldsymbol{x})=\Vert\boldsymbol{\hat{m}}(\boldsymbol{x})\Vert\Vert\boldsymbol{\hat{m}^{\dagger}}(\boldsymbol{x})\Vert, \end{equation}

where $\boldsymbol {\hat {m}}(\boldsymbol {x})$ and $\boldsymbol {\hat {m}^{\dagger} }(\boldsymbol {x})$ represent the direct and adjoint momentum vectors, respectively. In the previous section, the location of the wavemaker was identified by computing the intersection between the $\alpha ^+$ and $\alpha ^-$ branches (at $x/W=0.0731$). However, the exact location of the maximum structural sensitivity at this cross-flow plane is still unknown. This requires computation of the adjoint mode at this streamwise location. Therefore, the adjoint linear stability analysis is further pursued in this section.

Due to the inhomogeneous directions and hence differential dependencies in the linearized operator matrix, taking the Hermitian transpose of the direct operator as the adjoint operator, as has been done in Oberleithner, Rukes & Soria (Reference Oberleithner, Rukes and Soria2014); Müller et al. (Reference Müller, Lückoff, Paredes, Theofilis and Oberleithner2020), does not necessarily apply here. To find the adjoint of the system, the continuous approach is used. First, we define the inner product as

(4.20)\begin{equation} \left\langle\boldsymbol{\hat{q}}_1,\boldsymbol{\hat{q}}_2\right\rangle\equiv\int_{S}\boldsymbol{\hat{q}}_1^H\boldsymbol{\hat{q}}_2\,{{\rm d}}S\equiv\boldsymbol{\hat{q}}_1^H\boldsymbol{W}\boldsymbol{\hat{q}}_2. \end{equation}

The adjoint modes depend on this choice of norm, but when recombined with the direct modes to give the sensitivity measurement of the eigenvalue to changes in $\boldsymbol {\mathcal {L}}$, the effect of the norm cancels out (Chandler et al. Reference Chandler, Juniper, Nichols and Schmid2012).

By definition, the adjoint operator $\boldsymbol {\mathcal {L}^{\dagger} }$ should fulfil

(4.21)\begin{equation} \left\langle\boldsymbol{\hat{q}^{\dagger}},\boldsymbol{\mathcal{L}}\boldsymbol{\hat{q}}\right\rangle\equiv\left\langle\boldsymbol{\mathcal{L}^{\dagger}}\boldsymbol{\hat{q}^{\dagger}},\boldsymbol{\hat{q}}\right\rangle. \end{equation}

This definition is valid for any pair of vectors, but for convenience, they are expressed in terms of the direct state vector $\boldsymbol {\hat {q}}$ and the adjoint state vector $\boldsymbol {\hat {q}^{\dagger} }$, as done in Marquet et al. (Reference Marquet, Lombardi, Chomaz, Sipp and Jacquin2009); Qadri, Mistry & Juniper (Reference Qadri, Mistry and Juniper2013). More specifically, we consider the spatial stability form, then the generalized eigenvalue problem is pre-multiplied by the adjoint eigenvector

(4.22)\begin{equation} \left\langle\boldsymbol{\hat{q}^{\dagger}},\boldsymbol{\mathcal{A}}\boldsymbol{\hat{q}}\right\rangle-\left\langle\boldsymbol{\hat{q}^{\dagger}},\alpha\boldsymbol{\mathcal{B}}\boldsymbol{\hat{q}}\right\rangle=0. \end{equation}

By successively integrating the equation by parts using Green's theorem, (4.22) is rearranged to

(4.23)\begin{equation} \left\langle\boldsymbol{\mathcal{A}^{\dagger}}\boldsymbol{\hat{q}^{\dagger}},\boldsymbol{\hat{q}}\right\rangle-\left\langle\alpha^{\dagger}\boldsymbol{\mathcal{B}^{\dagger}}\boldsymbol{\hat{q}^{\dagger}},\boldsymbol{\hat{q}}\right\rangle=0. \end{equation}

Note that the integration-by-parts process would leave boundary terms in the expression. Therefore, appropriate adjoint boundary conditions are applied to cancel out all boundary terms. Then the adjoint of the direct eigenvalue problem is obtained from (4.23) as

(4.24)\begin{equation} \boldsymbol{\mathcal{A}^{\dagger}}\boldsymbol{\hat{q}^{\dagger}}=\alpha^{\dagger}\boldsymbol{\mathcal{B}^{\dagger}}\boldsymbol{\hat{q}^{\dagger}}. \end{equation}

The detailed derivations as well as the validation of the adjoint operators and boundary conditions used in the current study are presented in Appendix D.

It can be shown that, for well-converged adjoint solutions, the adjoint eigenvalues, ordered by the same rule as the direct ones, are the complex conjugates of the direct ones (Schmid & Henningson Reference Schmid and Henningson2000), and the bi-orthogonality condition writes

(4.25)\begin{equation} \left(\left(\alpha_m^{\dagger}\right)^\ast{-}\alpha_n\right)\left[\left(\boldsymbol{\hat{q}^{\dagger}}_m\right)^H\boldsymbol{W}\boldsymbol{\mathcal{B}}\boldsymbol{\hat{q}}_n\right]=0. \end{equation}

The bi-orthogonality enables the characterization of the receptivity of the open-loop forcing, and the sensitivity to modification of the mean flow (Marquet et al. Reference Marquet, Lombardi, Chomaz, Sipp and Jacquin2009).

In figure 16, the distribution of the structural sensitivity at the streamwise location where $\alpha ^+$ and $\alpha ^-$ branches intersect ($x/W=0.0731$) is shown. This field is obtained by first computing the adjoint mode on the considered cross-flow plane based on the methodology presented above and then calculating the $L^2$ norm of the sensitivity tensor following (4.19). In addition, mean flow streamlines are also included so that the structural sensitivity can be related to specific flow structures. Based on the results, both the upper and lower recirculation zones are highlighted. In particular, it can be observed that the position of the highest structural sensitivity, enclosed by the blue solid line in figure 16, is located slightly below the stagnation point of the recirculation region. The streamwise vortex pair generated from both sides of the train, on the other hand, seems to have no relevance to the sensitivity of the linear global mode. The highest structural sensitivity at this location, shown in figure 16, suggests that the interaction between the lower and upper shear layers, which drives the self-excitation of the instability, is the most sensitive to a steady external forcing and is therefore of significant importance in terms of controlling of the global mode. Note that, in this part, a deeper interpretation of different components of the structural sensitivity tensor, as has been done in Qadri et al. (Reference Qadri, Mistry and Juniper2013), is beyond the scope of this work.

Figure 16. Structural sensitivity indicating sensitivity to mean flow modifications; the values are normalized by the maximum value. Streamlines are drawn based on mean flow; the blue line indicates the normalized structural sensitivity = 0.90 isocontour.