2.1. Governing equations

The time-dependent three-dimensional flow field is obtained by solving the following governing equations. The continuity and momentum equations for a three-dimensional, incompressible and viscous flow are given as

(2.1)\begin{equation} \left.\begin{gathered} \boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{V} = 0, \\ \frac{\partial \boldsymbol{V}}{\partial t} =-\boldsymbol{\nabla} p -\frac{1}{2} \left[\boldsymbol{\nabla}(\boldsymbol{V}\otimes \boldsymbol{V}) + (\boldsymbol{V}\boldsymbol{\cdot}\boldsymbol{\nabla})\boldsymbol{V} \right]+\nu\, \nabla^2\boldsymbol{V} +f. \end{gathered}\right\} \end{equation}

In the above equation, $\boldsymbol {V}$ is the velocity vector with components $u$, $v$ and $w$ in the streamwise, wall-normal and spanwise directions, respectively. Here, $t$, $p$ and $\nu$ correspond to flow time, pressure and kinematic viscosity. The nonlinear terms in the governing equation are expressed in the skew-symmetric form since it is resilient to aliasing errors (Kravchenko & Moin Reference Kravchenko and Moin1997). The embedded body region in the computational domain is enforced by the body force field ($\,f$) in the momentum equation using the immersed boundary method (IBM). For numerically solving the governing equations, a high-order finite-difference flow solver Incompact3d (Laizet & Lamballais Reference Laizet and Lamballais2009; Laizet & Li Reference Laizet and Li2011), with a Cartesian mesh, is used. Incompact3d has been used extensively for various transitional and turbulent flow studies (Bempedelis & Steiros Reference Bempedelis and Steiros2022; Giri et al. Reference Giri, Biswas, Chase, Xue, Abkarian, Mendez, Saha and Stone2022).

In this code, spatial discretization of governing equations on a uniformly spaced Cartesian mesh is accomplished using a sixth-order compact finite-difference scheme. A third-order Adams–Bashforth approach is used for the time integration of the discretized governing equation. A staggered pressure grid from the velocity grid by half mesh length is implemented to avoid spurious pressure oscillations. The modified Poisson equation obtained by imposing the IBM is dealt with by spectral methods using the modified wavenumber formalism. 2Decomp&FFT (Li & Laizet Reference Li and Laizet2010), a domain decomposition library, performs fast Fourier transforms involved in spectral techniques. The library also contains a domain decomposition algorithm for efficient scaling and distribution of memory in high-performance computing systems.

2.2. Computational domain and boundary conditions

A sketch of the computational set-up of flow in a diverging channel is shown in figure 1. The computational domain chosen for this study is a small segment of the experimental set-up employed by Das et al. (Reference Das, Srinivasan and Arakeri2016). Here, $X$, $Y$ and $Z$ are streamwise, wall-normal and spanwise distances, respectively. In the simulation, the computational domain has length 1.2 m, width 0.142 m, which is equal to half the width of the experimental section, and height 0.15 m, as illustrated in figure 1(a). At the constant channel section, the embedded body section has height 0.07 m. After a length 0.3464 m, the edge starts to curve along an arc with radius (R) 0.1 m. Later, the curve joins smoothly to the diverging section with angle of depression $6.2 ^\circ$ (figure 1b). Similarly, the end of the diverging part joins fluidly with the bottom wall of the channel.

Figure 1. Computational domain along with boundary conditions: (a) three-dimensional view of the computational domain; (b) dimensions of the diverging section; (c) mean inflow variation during a pulse; and (d) variations in the inlet velocity profile during different velocity phases. (All the dimensions are in metres.)

A no-slip boundary condition is enforced on both top and bottom walls for devising identical experimental set-up conditions in the computational domain. A time-varying inlet condition based on the analytical solution of trapezoidal mean flow variation is imposed at the inlet of the computational domain. The free-slip condition is applied to both the right and left boundaries. The one-dimensional advective outflow equation is implemented as the exit boundary condition is given by

(2.2)\begin{equation} \left.\frac{\partial u}{\partial t} \right|_{{X}=l_x} + U_c \left.\frac{\partial u}{\partial {X}}\right|_{{X}=l_x} = 0, \end{equation}

where the mean velocity of the inlet velocity profile is taken as the advection velocity ($U_c$), and $l_x$ is the domain length in the streamwise direction.

The following equations give the mean inflow velocity for the four phases of piston motion:

(2.3)\begin{align} u_{p}(t)&=U_{p}\,\frac{t}{t_0 } \quad \text{for} \ 0\leq t\leq t_0, \nonumber\\ &=U_{p} \quad \text{for} \ t_0\leq t\leq t_1, \nonumber\\ &=U_{p}\,\frac{t_2-t}{t_2-t_1}\quad \text{for} \ t_1\leq t\leq t_2, \nonumber\\ &=0\quad \text{for} \ t> t_2. \end{align}

A trapezoidal pulse of mean inflow constitutes a constant acceleration phase ($0$ to $t_0$), a constant mean inflow phase ($t_0$ to $t_1$), and a constant deceleration phase ($t_1$ to $t_2$). Here, $U_p$ is the mean inlet velocity in the constant inflow phase (figure 1c). Such an inflow configuration can study the individual effects of the acceleration and deceleration phases of the inflow pulse. Analytical solutions of a trapezoidal mean inflow variation (2.3) for a two-dimensional channel following Das & Arakeri (Reference Das and Arakeri1998) are imposed at the inlet. Time-varying small-amplitude perturbations are generated by the combination of two parts of the analytical velocity solution (A5). The inlet boundary condition is not modified to account for free-stream turbulence. Details on the governing equations and solution procedures, along with velocity profiles imposed at the inlet during different phases, can be found in Appendix A. This approach shortens the entrance length and allows a shorter domain, reducing computational load. Since the inlet velocity profiles are built upon the assumption of a two-dimensional channel, using a subdomain guarantees the application of a slip boundary on the sidewalls. The work of Sarath & Manu (Reference Sarath and Manu2022) provides additional information about the procedure to obtain the inflow velocity profile.

Four typical velocity profiles imposed at the inlet during different phases are shown in figure 1(d). The analytical solution generates oscillating components, which cause small-amplitude oscillations, as shown in the figure. The amplitude of oscillations depends on the phase of the mean inflow velocity, and decreases when the mean inflow velocity is constant, as shown in instances V2 and V4. A spanwise width 0.142 m is sufficient to allow the domain to accommodate spanwise oscillations formed by the three-dimensional instability. Keeping the same channel height while reducing the area under consideration in the streamwise and spanwise directions allows for denser grid analysis and minimizes the computational load. A complete domain simulation demonstrated minor deviations (less than 1 %) in the wavelength of spanwise oscillation, affirming the selection of such a section for numerical analysis.

In the remaining sections of this paper, flow features are illustrated using non-dimensional spatial scales. Streamwise distance is non-dimensionalized as $x = ({X - X_s})/{h_b}$, where $h_b$ is the embedded body height at the inlet, and $X_s$ denotes the start of the diverging section ($X_s = 0.3464$). Wall-normal and spanwise distances are non-dimensionalized by using the body height defined by $y = {Y}/{h_b}$ and $z = {Z}/{h_b}$, respectively. As in the experiments of Das et al. (Reference Das, Srinivasan and Arakeri2016), the working fluid is selected to be water with kinematic viscosity ($\nu$) $10^{-6}\ \textrm {m}^2\ \textrm {s}^{-1}$. The parameters provided below are used to analyse the flow dynamics.

The Reynolds number is

(2.4)\begin{equation} Re_h = \frac{U_p h}{\nu}. \end{equation}

The acceleration Reynolds number is defined as

(2.5)\begin{equation} Re_a = \sqrt{\frac{a h^3}{\nu^2}}, \end{equation}

where $h$ is the inlet channel height, and $a$ is the acceleration (${U_p}/{t_0}$). Similarly, for varying deceleration cases, a deceleration Reynolds number is defined by

(2.6)\begin{equation} Re_d = \sqrt{\frac{d h^3}{\nu^2}}, \end{equation}

where $d$ is the deceleration (${U_p}/({t_2-t_1})$). In the present work, the Reynolds numbers based on the viscous length scales are defined as

(2.7a,b)\begin{equation} Re_{\delta} = \frac{U_p \delta}{\nu}\quad \text{and} \quad Re_{\delta^*} = \frac{U_p \delta^*}{\nu}, \end{equation}

where $\delta$ and $\delta ^*$ represents boundary layer and displacement thicknesses, respectively. Circulation of vortices for a particular time instance is calculated by $\varGamma _{\omega _z} = \iint _{A_{\varGamma }}\omega _z\,\textrm {d}A$. Here, the area $A_{\varGamma }$ is set appropriately to determine the circulation for top and bottom vortex flow features while omitting the wall boundary region (as depicted in figure 1b), and $\omega _z$ indicates the spanwise vorticity ($\omega _z = {\partial v}/{\partial X} - {\partial u}/{\partial Y}$). Similar to other vortex flow evolution studies (Le Dizes & Laporte Reference Le Dizes and Laporte2002; Leweke et al. Reference Leweke, Le Dizes and Williamson2016), the vortex Reynolds number based on spanwise circulation is estimated by

(2.8)\begin{equation} Re_{\varGamma_{\omega,z}} = \frac{\varGamma_{\omega_z}}{\nu}. \end{equation}

Table 1 shows the simulation parameters for 12 different flow cases. Each case is assigned an alphanumeric code to identify its simulation parameters. Of 12 simulations, cases with identical mean inflow velocity ($U_p$) are marked by letters A, B, C and D, respectively, for low, moderate, high and very high inflow velocities. In addition, numerals indicate cases with different deceleration parameters and the same Reynolds number: numbers 1, 2 and 3 refer to high, moderate and low deceleration. For all the cases, the acceleration Reynolds number ($Re_a$) is kept constant at 10 822. The Reynolds number based on the boundary layer thickness ($Re_{\delta _s}$) for all cases at separation time is given in table 1.

Table 1. Simulation parameters ($Re_a = 10\,822$).

As a measure of the spatial pressure distribution, the variation of the streamwise velocity component in the inviscid region ($U_x$) along the streamwise direction is shown in figure 2. The velocity profile is shown at half the constant velocity period. Since the acceleration and constant inflow velocity period remain the same, $\alpha$ and $\beta$ remain the same for varying deceleration cases. In figure 2, the continuous line represents the fitted curve, while the symbols represent the velocity obtained from the simulation. The symbols are placed at a distance of 15 grid point spaces between each pair. The velocity variation in the diverging section (from $X = 0.4$ to $X = 1.06$) can be approximated by using the relation

(2.9)\begin{equation} U_x = \alpha x^{\beta}, \end{equation}

where $\alpha$ and $\beta$ are constants and vary with Reynolds number ($Re_h$). The obtained values for $\alpha$ and $\beta$ for all cases are marked in figure 2.

Figure 2. Streamwise velocity variation in inviscid region along the streamwise direction ($y=1.4286$, $z=1.0$).

The following non-dimensionalized time scales are used to distinguish flow events. At first, flow time is non-dimensionalized by $t_2$ to differentiate between both the pulse phase and the zero mean inflow period ($t^* = {t}/{t_2}$). In order to compare different deceleration cases, a non-dimensionalized time scale is identified as $t_d^* = ({t-t_1})/({t_2-t_1})$. Also, a critical flow time associated with three-dimensionally unstable cases, $t^*_{3D}= ({t_{3D} - t_1})/({t_2-t_1})$, is identified from the $t_{3D}$ physical time at which a visible secondary instability initiates in three-dimensionally unstable cases. The $t^*_{3D}$ values observed for cases showing three-dimensional disintegration are provided in table 1. In low and moderate Reynolds number cases, the flow stays in the two-dimensional regime; hence $t^*_{3D}$ is absent for these cases.

2.3. Grid independence and numerical validation

The grid-independent analysis is performed by comparing the evolution of streamwise velocity (case D1) for different grids with elements $961 \times 193 \times 181$ (grid A), $1501 \times 301 \times 289$ (grid B) and $1921 \times 385 \times 361$ (grid C). Figure 3(a) shows the streamwise velocity through the central axis ($y=1.3143$, $z=1.0$) of the diverging channel at the end of the constant velocity phase for different grid sizes. The velocity component is non-dimensionalized by the maximum velocity magnitude for the flow instance ($u^* = {u}/{u_{max}}$). The velocity profiles for three different grid sizes are shown in figure 3(b) following the onset of initial instability at $x = 0.5$, $z = 1.0$. Grids A and B differ by approximately 1.5 %, three times greater than grids B and C. The root-mean-square (r.m.s.) deviation of streamwise velocity component ($u_{rms}$) of the velocity field during constant mean inflow phase is calculated as

(2.10)\begin{equation} u_{rms} = \left( \frac{1}{N} \sum_{l=1}^{l=N} (u'{(l)})^2 \right)^{1/2}, \end{equation}

where the velocity perturbation ($u'$) is calculated by

(2.11a,b)\begin{equation} u'(l) = u(l) - u_{mean}\quad \text{and}\quad u_{mean} = \frac{1}{N} \sum_{l=1}^{l=N} u{(l)}, \end{equation}

for $N$ snapshots belonging to the constant velocity phase. A comparison of the r.m.s. deviation of the streamwise velocity component developed over the constant mean inflow period reveals a difference of 8.4 % between grids A and B, while the difference between grids B and C is below 4 %. Consequently, grid B is selected for the numerical simulations due to accuracy and computational economy. Based on time step dependency analysis with time steps ranging from $1 \times 10^{-3}$ s to $1 \times 10^{-5}$ s, a time step $1 \times 10^{-4}$ s ($\textrm {CFL} = 0.02$) was found to be computationally and accurately affordable.

Figure 3. Computational model validation: (a) streamwise velocity variation in the streamwise direction (before separation, $t^* = 0.631$). (b) Streamwise velocity variation in the wall-normal direction (after separation, $t^* = 0.946$). (c) Experimental comparison of $t_s$.

Comparison of time of flow separation ($t_s$) with the results of Das et al. (Reference Das, Srinivasan and Arakeri2016) (figure 3b) validates the computational method. The previous experimental works incorporated a two-dimensional simulation of vortex evolution to calculate flow separation time. Here, by taking the flow time and position of zero wall shear stress, we determine the flow separation time and separation point. All the cases show an excellent match with the reported experimental values (less than 8 % difference). An experimental and simulation comparison of the flow field for a high-velocity case is provided in § 3.

5.1. The emergence of secondary instability and local breakdown

In high Reynolds numbers, the vortices formed by two-dimensional primary inflectional instability further undergo secondary instability, creating three-dimensional structures. Vortex flow structures developing in a high-deceleration case (C1) are revealed by the iso-surfaces of non-dimensionalized spanwise vorticity in figure 16. Primary inflectional instability causes the formation of negative vortices (LP1–LP4), which further induces secondary positive vortices (lp1–lp4) from the bottom wall boundary layer (figure 16a). This secondary vortex and the primary vortex form a pair near the wall proximity. Due to the mutual induction of the vortices, the pair detaches from the bottom wall, still pertaining to the two-dimensional nature ($t^* = 1.25$). Most flow features evolve near the initial diverging section emphasized by the non-dimensionalized streamwise scale on the top wall. The induced angular velocity of the upstream vortex pair (LP1, lp1) drives them to roll towards the vortex pair at the downstream location, (LP2, lp2). Such a flow development results in stretching and splitting of the positive secondary vortex by the primary vortices (lp2a, lp2b in figure 16b). The residual momentum pushes the secondary vortex (lp1) from the upstream pair to join the downstream couple; such a tri-vortex group further amplifies the roll-on process ($t^* = 1.93$). Inflectional profiles in the boundary layer seed multiple vortices from the boundary layer, as portrayed in figure 16(b) (pairs 5 and 6).

Figure 16. Temporal evolution of three-dimensional flow features identified by spanwise vorticity for case C1.

The secondary vortex structure (lp2b) exhibits a spanwise oscillation while orbiting the primary vortex (LP1), as shown in figure 16(c). Similar vortex flow features are exhibited by the vortex pairs detaching from the bottom wall surfaces downstream (pairs 3 and 4). A secondary vortex (lp4) downstream ($x\approx 1.0$) undergoes a similar fashion of vortex splitting, creating multiple positive vortices as in the former time instances for the upstream secondary vortex (lp2).

The merging of the primary vortex cores (LP4 and LP6) is visible in the same instance ($t^* = 1.93$). As the flow progresses, oscillations amplify in the secondary vortex circling the primary vortex. A sandwiching effect of the merged negative vortex cores (MG1 and MG2) stretches the secondary vortices around them (figure 16d). After the stretching, spanwise oscillations disintegrate into loops around the primary spanwise vortex flow structures (MG1 and MG2). However, the initially ejected secondary vortex (lp1) survives the vortex interactions and shows a three-dimensional oscillation while orbiting the primary vortex. In a later flow instance ($t^* = 2.72$), the secondary vortex structure breaks down into spanwise loops around partially disintegrating primary vortices (LP1; figure 16e). A complete transition to a turbulent structure with a negative vorticity core is observed as the flow progresses ($t^* = 3.18$). Simultaneous production of two-dimensional vortices and merging transitions is observable in the upstream and downstream positions in figure 16(f). Since this flow development lies in the zero mean inflow phase, advection of three-dimensional roll-ups is not evidenced for case C1, and flow features are confined to a small area near the initial diverging section ($x = 0.4- 1.5$). For case C1, flow evolution over the top wall has not yet been developed at $t^* = 3.18$. Induction from the bottom wall vortices creates a positive vorticity patch over the top wall at $t^* = 3.18$, indicating vortex roll-up in later flow time.

The three-dimensional topology of vortex structures is identified further using the $\lambda _2$ method and contoured by non-dimensional streamwise vorticity in figure 17, which illustrates the spanwise oscillation developed over the secondary vortex for case C1 at $t^* = 1.93$. The alternate streamwise vorticity originated over the spanwise vortex roll-ups displays a symmetrical pattern. A cross-section ($Y$–$Z$) through the central plane of vortex pairs ($x = 0.308$) is taken, to understand the nature of oscillations. Non-dimensionalized spanwise vorticity contoured snapshots (figure 17b) provide a clear visual of secondary vortex developing oscillations. The presence of alternating streamwise vorticity is further explained by superimposing the vector plot in figure 17(c). Circulation regions are discernible over the spanwise vortex roll-ups, as seen in the inset.

Figure 17. (a) Flow features identified (LP1 and lp2b) using the $\lambda _2$ method contoured by streamwise vorticity ($\lambda _2 = -1$). (b) Section through the spanwise oscillation contoured by spanwise vorticity ($x = 0.29$). (c) Section through the spanwise oscillation contoured by streamwise vorticity.

The temporal evolution of flow vortices in cases similar to case C1 (type II) is evidenced in figure 18. Vortex flow structures at two flow instances are shown: the first flow instance portrays the three-dimensional oscillations during the early transition phase, while the second instance illustrates the flow structures in the turbulent phase. Identical to case C1, primary instability in the deceleration phase develops into a vortex pair, which moves downstream and interacts with similar vortex pairs. Figures 18(a,b) indicate the formation of spanwise oscillations over the secondary vortices in case C2. Stretching of the secondary vortex around the primary vortex is evidenced at $t^* = 2.0$. Along with three-dimensional disintegration, vortex pairs eject into the flow core at a downstream position due to the induced velocity. As flow progresses, analogous to case C1, the formation of a turbulent three-dimensional structure with negative spanwise vorticity is illustrated in figures 18(a,b). Unlike case C1, the formation of top wall vortices is evidenced for case C2 as depicted in figures 18(a,b). During the zero mean inflow phase at $t^* = 2.13$, the ejected vortex pair interacts with the top wall boundary layer. Figures 18(c,d) and 18(e,f) show identical flow evolution for cases D1 and D2. Three-dimensional disintegration initiates with vortex pair interaction, resulting in a locally turbulent region. The turbulent region develops earlier for high Reynolds number cases (D1 and D2) than for high flow Reynolds number cases (C1 and C2).

Figure 18. Temporal evolution of three-dimensional flow features identified by spanwise vorticity for cases belonging to type II instability: (a,b) case C2, (c,d) case D1, and (e,f) case D2.

5.2. Modelling of secondary instability

The DMD algorithm is used to analyse secondary instability over the shear layer vortices resulting from primary inflectional instability; it identifies the underlying dynamics of the coherent flow features developed during the flow evolution. We use the snapshot-based approach introduced by Schmid (Reference Schmid2010) to identify the secondary instability features along with their temporal dynamics. A review by Taira et al. (Reference Taira, Hemati, Brunton, Sun, Duraisamy, Bagheri, Dawson and Yeh2020) offers further examples of this approach to identify underlying features in many transitional and turbulent flows.

The DMD algorithm starts by arranging state vectors from simulation data snapshots column-wise into two snapshot matrices $\boldsymbol {U}_1$ and $\boldsymbol {U}_2$, with two consecutive time intervals ($\Delta \tau$). DMD analysis identifies a best-fit linear operator $\boldsymbol {A}$ that relates matrix $\boldsymbol {U}_1$ to the matrix $\boldsymbol {U}_2$:

(5.1)\begin{equation} \boldsymbol{U}_2 = \boldsymbol{A} \boldsymbol{U}_1. \end{equation}

By taking the singular value decomposition (SVD) of the snapshot matrix ($\boldsymbol {U}_1$), we obtain left and right eigenvectors ($\boldsymbol {S}$ and $\boldsymbol {D}$, respectively), along with the eigenvalue ($\varLambda$). In order to build a best-fit linear operator correlating both matrices, we take the pseudo-inverse of $\boldsymbol {U}_1$ (by taking the conjugate transposes of SVD vectors ($\boldsymbol {S}^*,\boldsymbol {D}^*$) together with $\boldsymbol {U}_2$ as follows:

(5.2)\begin{equation} \boldsymbol{\mathsf{A}} = \boldsymbol{U}_2 \boldsymbol{D}^* \varLambda^{-1} \boldsymbol{S}^*. \end{equation}

By eigen-decomposition, the low-dimensional model is decomposed into eigenvector ($\boldsymbol {W}$) and eigenvalues ($\boldsymbol {\lambda }_D$). The temporal dynamics of the system can be identified from the eigenvalues using the following relations for growth rate and frequency:

(5.3)$$\begin{gather} \sigma_{DMD} = \frac{\log\left({\rm Re}(\lambda_D)\right)}{2 {\rm \pi}\, \Delta \tau}, \end{gather}$$

(5.4)$$\begin{gather}f_{DMD} = \frac{\log\left({\rm Im}(\lambda_D)\right)}{2 {\rm \pi}\,\Delta \tau}. \end{gather}$$

In order to obtain the DMD modes ($\varPhi$), the low-dimensional model is reconstructed from its eigenvectors as

(5.5)\begin{equation} \boldsymbol{\varPhi} = \boldsymbol{U}_2 \boldsymbol{D}\varLambda^{-1}\boldsymbol{W}. \end{equation}

Through a three-dimensional DMD of streamwise vorticity data, the most unstable flow features associated with secondary instabilities are identified (figure 19). For the DMD calculation, three-dimensional snapshots are taken with time step $\Delta \tau =0.05$ s between each snapshot. The Ritz circle obtained from the DMD for cases belonging to the second category is presented in figure 19(a). The position of the mode with respect to the unit circle outlined in the figure indicates the stability of the modes. In general, a mode lying outside a circle indicates an unstable mode, while lying within signifies a stable mode; and when it lies on a circle, it is neutrally stable. In all cases, at least one mode displays an unstable trait, illustrated by the circle in figure 19(a). Figure 19(b) shows the growth rate and frequency distribution for the DMD modes. The highest growth rate mode obtained from DMD analysis is indicated by $1_{\sigma }$.

Figure 19. Three-dimensional DMD analysis results: (a) Ritz circle, (b) growth rate versus frequency, and (c) leading secondary instability modes based on growth rate criterion for type II cases.

When streamwise vortex growth is unstable, vortex roll-ups will experience three-dimensional destabilization, as manifested by the positive growth rates for all cases. Unstable three-dimensional modes are identified through the growth rate criteria for all cases and are shown in figure 19. The three-dimensional morphology of unstable modes is visualized using an iso-surface of streamwise vorticity. A spanwise variation of the streamwise vorticity is evidenced in all cases. The spanwise wavelength obtained for the mode with the highest growth rate is provided in table 2. Averaging the distance between the peaks of the streamwise vorticity plotted across the vortex core provides the mean spanwise wavelength. The average spanwise wavelength observed for all three-dimensional cases is tabulated in table 2. Both wavelengths are non-dimensionalized by the distance between the cores of vortex pairs ($b$ is the distance between maximum and minimum vorticity magnitudes). The spanwise wavelength for the coherent flow features identified by DMD analysis ($\lambda _{DMD}$) lies close to the mean spanwise wavelength determined from streamwise vorticity variation over the oscillation ($\lambda _{mean}$).

Table 2. Spanwise wavelength comparison and parameters for Lamb–Oseen approximation.

The secondary instabilities are formed during the zero mean phase in cases belonging to the type II category, which is similar to the short-wavelength elliptic instability demonstrated by Laporte & Corjon (Reference Laporte and Corjon2000). Since multiple vortex pairs are observed near the bottom wall, unlike classical short wavelength vortex instability, theoretical stability analysis is performed for most magnified vortex pairs. Variation of the spanwise vorticity in a vortex pair ejecting from the bottom wall (case C1) is shown in figure 20(a). An identical depiction of an approximated Lamb–Oseen pair obtained from (5.6) is given in figure 20(b). The vorticity distribution for such a vortex pair with circulations $\varGamma _1$ and $\varGamma _2$ may be approximated using the Lamb–Oseen equation (Leweke et al. Reference Leweke, Le Dizes and Williamson2016)

(5.6)\begin{equation} \omega_z = \frac{\varGamma_1}{{\rm \pi} a_1^2}\exp\left(-\frac{r^2}{a_1^2}\right) + \frac{\varGamma_2}{{\rm \pi} a_2^2}\exp\left(-\frac{r^2}{a_2^2}\right), \end{equation}

where $r$ represents the distance from the vortex core position. The core radius ($a_1$ and $a_2$) for a vortex centred at $\boldsymbol {X}^c$ is obtained from

(5.7)$$\begin{gather} (x_{c_1},y_{c_1}) = \left(\frac{1}{\varGamma_1}\int_{D_1} x \omega_z\, {\rm d}S, \ \frac{1}{\varGamma_1}\int_{D_1} y \omega_z\,{\rm d}S\right), \end{gather}$$

(5.8)$$\begin{gather}a_1^2 = \frac{1}{\varGamma}\int_{D_1}|\boldsymbol{X}-\boldsymbol{X}^c|^2 \omega_z \,{\rm d}S, \end{gather}$$

where $D_1$ represents the domain containing each vortex, $\omega _z$ defines the spanwise vorticity, and $\textrm {d}S$ represents the infinitesimal area.

Figure 20. Depiction of Lamb–Oseen model: (a) simulation, (b) Lamb–Oseen model, and (c) Lamb–Oseen model comparison for vortex pair in locally unstable cases.

To ascertain the validity of the Lamb–Oseen approximation, the plot of the spanwise vorticity obtained from (5.6) is compared with the simulation results (through the vorticity cores). Table 2 provides the essential parameters obtained from simulation used to estimate the Lamb–Oseen approximation for locally unstable cases (type II) as depicted in figure 20(c). Here, $b_{LO}$ denotes the distance between the first positive vortex core and the second negative vortex in approximated vorticity distribution. In cases C1, C2 and D2, a vertical pair of vortices is compared with a vertical Lamb–Oseen pair, while in case D1, a better approximation is obtained for a horizontal pair and is compared to a similar Lamb–Oseen approximation (figure 20c). As presented in figure 20, the assumption of a Lamb–Oseen vortex approximation remains true in all selected flow instances.

Le Dizes & Laporte (Reference Le Dizes and Laporte2002) proposed an explicit relation using approximate linear expressions for the internal strain ratio ($s_r$) and inertial wave vector inclination for predicting the growth rate of an elliptic instability in a counter-rotating vortex pair. Leweke et al. (Reference Leweke, Le Dizes and Williamson2016), in their review, presented a revised linear fit to determine the frequencies ($\omega$) and damping rates ($\zeta$) of the first two Kelvin modes. The growth rate for the first two modes of elliptic instability in a Lamb–Oseen vortex pair is given by

(5.9)\begin{align} \sigma_1^{*,{(n)}} = \sqrt{\left(\frac{3}{4} - \frac{\bar{\varOmega}_1}{4}\right)^4 s_r^2\left(\bar{\varOmega}_1\right)\left(\frac{\varGamma_2}{\varGamma_1}\right)^2 - \left(\omega^{(n)}- \bar{\varOmega}_1\right)^2 \left(\frac{b}{a_1}\right)^4}- \left(\frac{b}{ a_1}\right)^2 \left(\frac{\zeta^{(n)}}{Re_{\varGamma_1}}\right) . \end{align}

The superscript $n$ represents the mode number, and the subscript denotes the vortex number. Here, the growth rate of the mode is non-dimensionalized by the time scale of translational motion (${2{\rm \pi} b^2}/{\varGamma _1}$). Linear expressions for real ($\omega$) and imaginary ($\zeta$) parts of the complex frequency in (5.9) are

(5.10a)$$\begin{gather} \bar{\varOmega}_1 = \left(\frac{a_1}{b}\right)^2\left(\frac{\varGamma_1 + \varGamma_2}{\varGamma_1}\right) , \end{gather}$$

(5.10b)$$\begin{gather}s_r(\bar{\varOmega}) = 1.5+0.1323\left(0.32-\bar{\varOmega}\right)^{-9/5} , \end{gather}$$

(5.10c)$$\begin{gather}\omega^{(1)} =-0.135\left(ka_1 - 2.26\right), \end{gather}$$

(5.10d)$$\begin{gather}\omega^{(2)} =-0.084\left(ka_1 - 3.95\right), \end{gather}$$

(5.10e)$$\begin{gather}\zeta^{(1)} = 74.02 + 64.15\left(ka_1 - 2.26\right) , \end{gather}$$

(5.10f)$$\begin{gather}\zeta^{(2)} = 229.6 + 104.3\left(ka_1 - 3.95\right) . \end{gather}$$

Using (5.9), the growth rate of elliptical instability for a counter-rotating vortex pair ejected from the bottom boundary of the wall was determined. The growth rate curves obtained for different spanwise wavenumbers for secondary vortices are presented in figure 21. A similar analysis was conducted on the vortex structures in type I cases. The secondary vortices, which develop in two-dimensional advecting and diffusing cases, indicate stable growth rates and have been confirmed by our calculations but are not shown. The growth rate curves are calculated for the secondary vortex, which undergoes three-dimensional disintegration (using Lamb–Oseen approximation as shown in figure 20c). The growth rate curves indicate an unstable first mode alongside a stable second mode (figure 21a). An identical growth rate curve is obtained for case D1 (figure 21b), revealing the formation of first-mode elliptic instability in high-deceleration cases (C1 and D1). Unlike high-deceleration cases, both the first and second modes show an unstable nature in moderate-deceleration cases (C2 and D2) (figures 21c,d).

Figure 21. Growth rate curves for first and second elliptic modes for secondary vortex in locally unstable cases: (a) case C1, (b) case D1, (c) case C2, and (d) case D2.

The temporal growth of vorticity/spanwise oscillation can be associated with the secondary instability (elliptic instability) of the shear layer vortices formed by the primary instability mechanism. The small spanwise velocity component can grow into a larger magnitude if the flow conditions are susceptible to elliptical instability. The fluctuations amplify due to the elliptic instability to a significant value ($10^{-2}$) for all type II cases, as is evident in figure 15, whereas it is negligible for type I cases. We also found that the amplitude of spanwise fluctuation is higher for case D2, showing a maximum growth rate. Identical to the growth rate curves, case D2 indicates a maximum amplitude of order 0.01. Compared to very high Reynolds cases (D1 and D2), high Reynolds cases (C1 and C2) possess lower amplitude, and the growth rate for case C1 shows the lowest value.

6.1. Unsteady separation and flow breakdown

Vortex flow structures identified by the non-dimensional spanwise vorticity in case C3 are portrayed in figure 23. Shear layer roll-ups marked by BS1–BS5 are observed at nearly equally spaced locations over the diverging section. Vortex flow structures are highlighted according to their formation sequence during flow evolution. Due to unsteady flow separation, vortices develop over the initial diverging section (PS1) during deceleration. Vortices evolving over the diverging section advect downstream, creating additional positive vortices (bs1) from the wall surface ($t^* = 0.55$). Persistent streamwise velocity deceleration leads to a continuous shedding of vortices from the initial diverging section (PS1 and PS2). A higher advective velocity of the upstream vortex leads to the pairing of the vortex roll-ups (BS1 and BS3, $t^* = 0.63$). Due to mutual induction, positive vortices (ps1a–ps1c) eject from the bottom wall vorticity layer and revolve around the primary vortex (PS1) while advecting downstream ($t^* = 0.69$). In such advecting vortex pairs, secondary vortices generate spanwise oscillations (ps1c) similar to locally unstable cases (type II). The formation of three-dimensional oscillations is attributed to the local interaction of advecting pairs. Simultaneous merging of vortex flow features (BS2 and BS4) and the three-dimensional disintegration of secondary vortices (ps1c, ps2a and ps2b) are evident in later deceleration flow instances ($t^* = 0.74$ and $0.76$).

Figure 23. Temporal evolution of three-dimensional flow features identified by non-dimensional spanwise vorticity for case C3.

A complex structure with multiple small-scale vortices is formed due to the growth of the three-dimensional oscillations into tube-like structures (MG4), which displace downstream due to the streamwise velocity ($t^* = 0.81$). Three-dimensional oscillations develop over the unsteady separation vortices when the flow passes half of the deceleration phase. A vortex pair discharging from the initial diverging section displays three-dimensional oscillation after the ejection ($t^* = 0.91$). Advecting vortex roll-ups interact with other flow features, generating a turbulent flow structure ($t^* = 0.91$). The decay of streamwise velocity further results in flow detachment, leading to a turbulent flow evolution during the zero mean inflow phase. Analogous to the former category, the flow features move backwards as the flow progresses in the zero mean inflow region. Tiny loop structures are formed through the disintegration of vortex flow structures at later flow instances ($t^* = 1.18$).

Identical to the bottom wall, the top wall boundary layer displays vortex roll-ups and three-dimensional disintegration, as presented in figure 23. An extended stay of the bottom wall vortex in the initial diverging section induces the flow to separate from the top wall leading to vortex roll-up. The top wall boundary layer broadens over the bottom separation region once the flow separates over the bottom wall ($t^* = 0.55$). In type III cases, flow separation occurs earlier in the deceleration phase, forcing the fluid over the top wall. Formation of a primary positive top wall vortex (TS1) subsequently results in the production of a negative secondary vortex (ts1, $t^* = 0.69$). The shedding of vortices due to unsteady separation promotes the development of flow structures (TS pairs 2 and 3) across the top wall ($t^* = 0.74$). The top wall vortices retain their two-dimensional traits, while the bottom wall vortex pair undergoes three-dimensional disintegration at $t^* = 0.76$. Due to induced rotation, primary and secondary top wall vortices eject from the top wall boundary layer ($t^* = 0.81$). Analogous to the bottom wall flow features, top wall vortices clearly show the pairing behaviour at $t^* = 0.91$ (TS3). Towards the end of the deceleration phase, the top wall formations also generate three-dimensional oscillations. In the initial zero mean inflow phase, the interaction between the top and bottom wall structures is minimal ($t^* = 1.02$), while at a later flow instance ($t^* = 1.18$), the mixing of flow features over the top and bottom walls results in a turbulent flow evolution.

Similar to case C3, vortex evolution in case D3 demonstrates vortex generation and three-dimensional breakdown of flow characteristics during the deceleration phase (figure 24). The inflectional instability in the boundary layer develops into vortex roll-ups (BS1–BS5) during the deceleration phase ($t^* = 0.43$ and $0.49$). Compared to case C3, a higher mean inflow velocity results in a higher advective velocity for the flow features along with the vortex shedding (PS pairs) due to unsteady separation ($t^* = 0.52$). Three-dimensional oscillations induced over the secondary vortices amplify and disintegrate at subsequent flow instances as vortex pairs are ejected from the initial diverging section. In case D3, the vortex flow patterns disintegrate at an early pulse stage, resulting in an advecting turbulent structure that stays closer to the bottom wall ($t^* = 0.6, 0.64$). Unlike case C3, the top wall vortices (TS pairs 1 and 2) advect downstream along with the bottom wall structures due to the earlier inception of top wall structures. Flow formations over the top wall retain their two-dimensional nature at this flow instance ($t^* = 0.64$). As the flow moves forward through the deceleration phase, the advecting top wall vortices (TS1, TS2) evidence three-dimensional disintegration ($t^* = 0.7$). The turbulent flow features over the top and bottom walls advect downstream as flow forwards, while the interaction of the structures results in a turbulent flow as evidenced at $t^* = 0.75$ and $0.82$.

Figure 24. Temporal evolution of three-dimensional flow features identified by non-dimensional spanwise vorticity for case D3.

A separation bubble forms near the initial diverging region due to unsteady separation, constituting a spanwise vortex roll and induced positive vortex in the boundary layer. A close-up image of the vortex shedding and three-dimensional disintegration of the separation bubble vortices is illustrated in figure 25. Similar to the shedding process observed by Wissink & Rodi (Reference Wissink and Rodi2006), consecutive formation of the vortex rolls is evidenced during the deceleration phase (figure 24). As the flow decelerates further, the streamwise velocity weakens, amplifying the perturbations developing due to vortex interactions. As the flow decelerates ($t^* = 0.61$), the secondary vortices generate a spanwise oscillation during the shedding process, and develop into spanwise loops over the primary negative vortex (PS2). Decay in streamwise velocity induces an oscillation in the separation bubble due to the amplification of perturbations (PS3). As the flow progress ($t^* = 0.63$), the vortex structure (PS3) sheds downstream, creating a short interval and a turbulent separation bubble (PS4). Further deceleration leads to the disintegration of the separation bubble and moves to a turbulent regime ($t^* = 0.65$). The formation of a turbulent separation bubble leads to an uneven shedding of the turbulent flow features similar to the flow evolution observed by Wissink & Rodi (Reference Wissink and Rodi2006) ($t^* = 0.68$).

Figure 25. Three-dimensional disintegration of the separation bubble.

6.2. Shear layer shedding characteristics

In figure 26, the spanwise vorticity at a downstream location is probed to identify the temporal characteristics of the vortex shedding due to unsteady separation for cases C3 and D3. In case C3, vortex generation initiates around $t^* = 0.3$, which is at an earlier flow instance than in case D3. Shortly after the vortex ejection, the formation of a positive vortex is evident by a sharp positive peak. Further oscillations are indicative of the subsequent formation and shedding of vortex structures. High-frequency perturbations in the zero mean inflow phase are triggered by three-dimensional fragmentation of the flow features in case C3. The vortex shedding in case D3 is more frequent and reveals a three-dimensional breakdown at an earlier stage in the flow. The underlying frequency of the vortex shedding is identified by the frequency spectra obtained from the Fourier transforms of the vorticity evolution. The highest peak in the frequency spectra obtained by the Fourier analysis indicates the shedding frequency and is marked in figure 26(b) for cases C3 and D3. A second dominant frequency lies close to the subsequent harmonics of the preceding dominant frequency in both cases.

Figure 26. Spectra analysis of vorticity probe for spatially unstable cases. (a) Temporal variation of spanwise vorticity (case C3, $x = 2.633$, $y = 0.914$; case D3, $x = 2.32$, $y = 0.964$). (b) Frequency spectra of the spanwise vorticity variation.

The temporal characteristics of the periodic shedding of two-dimensional flow structures in type III cases are revealed by the DMD analysis of two-dimensional snapshots of spanwise vorticity. Unlike the former category cases, rapid evolution and the streamwise advection of flow features due to the low deceleration rate hinders the three-dimensional DMD analysis of streamwise vorticity. A total of 75 snapshots lying in the initial deceleration period between $t_i$ and $t_f$ are used for the DMD analysis for both cases, with time step $\Delta \tau =2 \times 10^{-2}$. Frequency was invariant when the number of snapshot sizes was increased. Figure 27 summarizes the DMD analysis of the spanwise vorticity in spatially unstable cases. Typically, the Ritz circles plotted in figures 27(a) and 27(b) indicate the stability of DMD modes. The red symbols indicate the growing modes of cases C3 and D3 that lie outside the unit circle. To further emphasize the destabilizing nature, figures 27(c) and 27(d) present bar diagrams of frequency ($\,f_{DMD}$) plotted against the growth rate ($\sigma$) determined using DMD analysis. In each case, $1_{\sigma }$ denotes the most unstable mode based on the growth rate criteria. Growth rates and frequency information for each mode are provided in each figure. Figures 27(e) and 27(f) display the most prominent (highest growth rate) mode for cases C3 and D3, respectively. The peak frequencies of both modes are higher in the high Reynolds number case D3. Invariably, the first mode represents the growth of the vortex from the separation bubble over the diverging channel. A second dominant mode holds a frequency nearly equal to the second harmonic, indicating the second dominant peak of the vorticity frequency spectra (figures 27g,h). An alternate pattern in the second mode indicates the spanwise vortex roll-ups over the diverging region. In case D3, early vortex development leads to continuous shedding mode, as shown in figure 27(h). DMD results also point to the vortex shedding from the separation bubble in this period. Both analyses indicate two-dimensional vortex shedding characteristics, which later disintegrate three-dimensionally during the zero mean inflow phase.

Figure 27. DMD analysis of spanwise vorticity evolution: (a) Ritz circle (case C3), (b) Ritz circle (case D3). Frequency versus growth rate: (c) case C3, (d) case D3. (e) First mode, case C3, (f) first mode, case D3, (g) second mode, case C3, and (h) second mode, case D3).

A consolidated comparison of temporal characteristics obtained from frequency spectra with the results from DMD analysis is included in table 3. An average displacement thickness is calculated for the velocity profiles over the separation point and is denoted as $\delta _{avg}^*$. Similarly, an average velocity is calculated by taking the mean of the average velocity of velocity profiles over the separation point (between $t_i$ and $t_f$). The Strouhal number ($St = {f\delta _{avg}^*}/{U_{avg}}$) is calculated for both probe analysis and DMD analysis by using the average displacement thickness along with the average velocity. In the present study, the unstable mode frequency scales with the viscosity length scales generally used in boundary layer transition studies (Klebanoff, Cleveland & Tidstrom Reference Klebanoff, Cleveland and Tidstrom1992; Bakchinov et al. Reference Bakchinov, Grek, Klingmann and Kozlov1995), and identified Strouhal frequencies lie near $0.2$. In light of the fact that only 2 out of 12 cases show periodic vortex shedding, the results are insufficient to support the generalization of the Strouhal number relation.

Table 3. Temporal characteristics of vortex shedding due to unsteady separation.