2.1 Problem formulation

Flow evolution is described by the Navier–Stokes and continuity equation of the form

(2.2) $$\begin{eqnarray}\displaystyle \left.\begin{array}{@{}c@{}}\displaystyle \unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}=0,\\ \displaystyle \frac{\unicode[STIX]{x2202}\boldsymbol{u}}{\unicode[STIX]{x2202}t}+\boldsymbol{u}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\boldsymbol{u}=-\unicode[STIX]{x1D735}p+\frac{1}{Re}\unicode[STIX]{x1D6FB}^{2}\boldsymbol{u},\end{array}\right\} & & \displaystyle\end{eqnarray}$$

where $\boldsymbol{u}=[u,v,w]^{\text{T}}$ is the velocity vector, $p$ is the pressure, the maximum of the streamwise velocity in the reference smooth channel, $U_{s}$ , has been used as the velocity scale, $\unicode[STIX]{x1D71A}U_{s}^{2}$ has been used as the pressure scale with $\unicode[STIX]{x1D71A}$ denoting density, $h/U_{s}$ plays the role of the time scale and $Re=(U_{s}h)/\unicode[STIX]{x1D708}$ stands for the Reynolds number with $\unicode[STIX]{x1D708}$ being the kinematic viscosity. The boundary conditions at the walls have the form of

(2.3) $$\begin{eqnarray}\displaystyle \boldsymbol{u}=0\quad \text{at }\left\{\begin{array}{@{}l@{}}y=y_{u},\quad \\ y=y_{l}.\quad \end{array}\right. & & \displaystyle\end{eqnarray}$$

We shall impose periodicity conditions in the $x$ -direction as a periodic geometry is expected to create periodic flow modulations. In general, periodicities of geometry and flow modulations do not need to be identical and this fact needs to be accounted for while selecting the computational domain. Since this domain needs to be finite, direct analysis of incommensurable states is not possible but the potential existence of such states could be inferred by using limiting processes and analytical arguments.

The reference flow has the form of plane Poiseuille flow confined between flat rigid walls at $y=\pm 1$ . The fluid motion has velocity $\boldsymbol{u}$ , pressure $p$ and flow rate $Q_{r}$ of the form:

(2.4) $$\begin{eqnarray}\displaystyle \left.\begin{array}{@{}c@{}}\displaystyle \boldsymbol{u}=[1-y^{2},0],\\ \displaystyle p=-\frac{2x}{Re},\\ \displaystyle Q_{r}=\frac{4}{3}.\end{array}\right\} & & \displaystyle\end{eqnarray}$$

It is assumed that while the channel morphs into the converging–diverging form, the pressure gradient driving the flow remains the same, resulting in changes of the flow rate. This process corresponds to the imposition of the pressure-gradient constraint of the form

(2.5) $$\begin{eqnarray}\displaystyle \left.\frac{\unicode[STIX]{x2202}p}{\unicode[STIX]{x2202}x}\right|_{mean}=\frac{-2}{Re}, & & \displaystyle\end{eqnarray}$$

with the resulting flow rate determined a posteriori. This means that the effective Reynolds number is a function of the groove geometry.

Problem (2.2)–(2.5) is solved numerically using an incompressible Navier–Stokes solver available in the NEKTAR $++$ software package (Cantwell et al. Reference Cantwell, Moxey, Comerford, Bolis, Rocco, Mengaldo, Grazia, Yakovlev, Lombard and Ekelschot2015; Bolis et al. Reference Bolis, Cantwell, Kirby and Sherwin2014). Spatial discretization is based on the spectral element method. Simulations have been carried out using a regular, structured quadrilateral mesh combined with fifth-order polynomials. The mesh has been generated using the GMSH package (Geuzaine & Remacle Reference Geuzaine and Remacle2009). The computational box extends over an integer number of corrugation wavelengths in the streamwise direction. Selection of the number of wavelengths is dictated by the flow physics. The flow is assumed to be periodic in the streamwise and spanwise directions. The temporal discretization uses the second-order velocity-correction scheme (Guermond & Shen Reference Guermond and Shen2003). The kinetic energy evolution criterion is used to test the existence of steady states. Average kinetic energy per unit area over the computational domain $\unicode[STIX]{x1D6FA}$ is defined as

(2.6) $$\begin{eqnarray}\displaystyle E_{k}(t)=\frac{1}{\unicode[STIX]{x1D6FA}}\int _{\unicode[STIX]{x1D6FA}}\frac{1}{2}\Vert \boldsymbol{u}(t)\Vert ^{2}\,\text{d}\unicode[STIX]{x1D6FA} & & \displaystyle\end{eqnarray}$$

and the simulations are assumed to reach a steady state when the energy change $\unicode[STIX]{x1D6FF}E_{r}$ defined as:

(2.7) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6FF}E_{r}=\frac{\underset{T_{f}-T_{l}\leqslant t\leqslant T_{f}}{\max }E_{k}(t)-\underset{T_{f}-T_{l}\leqslant t\leqslant T_{f}}{\min }E_{k}(t)}{\underset{T_{f}-T_{l}\leqslant t\leqslant T_{f}}{\max }E_{k}(t)} & & \displaystyle\end{eqnarray}$$

reaches a certain, predefined threshold in a suitably selected time interval. Here, $T_{f}$ represents the simulation final time and $T_{l}$ stands for the time interval over which the stationarity test is conducted. In the current study we use $T_{l}=10^{4}$ and $\unicode[STIX]{x1D6FF}E_{r}=10^{-6}$ . Verification of the spatial accuracy and testing required for the grid selection have been carried out using a computational box containing one corrugation wavelength with details given in appendix A. If the computational box contains several corrugation wavelengths, the same grid is used in each section corresponding to one corrugation wavelength.

2.2 Description of two-dimensional states

2.2.1 Long-wavelength corrugation

The stationary state in the limit of long corrugation wavelength ( $\unicode[STIX]{x1D6FC}\rightarrow 0$ ) can be determined analytically, following Mohammadi & Floryan (Reference Mohammadi and Floryan2012). The first step is to regularize the flow domain using a transformation of the form:

(2.8) $$\begin{eqnarray}\displaystyle \left.\begin{array}{@{}c@{}}\displaystyle \unicode[STIX]{x1D709}=\unicode[STIX]{x1D6FC}x,\\ \displaystyle \unicode[STIX]{x1D702}=\frac{y-1+S\cos (\unicode[STIX]{x1D6FC}x)}{1-S\cos (\unicode[STIX]{x1D6FC}x)}+1\end{array}\right\} & & \displaystyle\end{eqnarray}$$

leading to the field equations of the form:

(2.9) $$\begin{eqnarray}\displaystyle \left.\begin{array}{@{}c@{}}\displaystyle \unicode[STIX]{x1D6FC}\frac{\unicode[STIX]{x2202}u}{\unicode[STIX]{x2202}\unicode[STIX]{x1D709}}+F_{11}\frac{\unicode[STIX]{x2202}u}{\unicode[STIX]{x2202}\unicode[STIX]{x1D702}}+F_{12}\frac{\unicode[STIX]{x2202}v}{\unicode[STIX]{x2202}\unicode[STIX]{x1D702}}=0,\\ \displaystyle -\frac{\unicode[STIX]{x2202}^{2}u}{\unicode[STIX]{x2202}\unicode[STIX]{x1D702}^{2}}-F_{1}\frac{\unicode[STIX]{x2202}u}{\unicode[STIX]{x2202}\unicode[STIX]{x1D702}}+F_{2}u\frac{\unicode[STIX]{x2202}u}{\unicode[STIX]{x2202}\unicode[STIX]{x1D702}}+F_{3}v\frac{\unicode[STIX]{x2202}u}{\unicode[STIX]{x2202}\unicode[STIX]{x1D702}}-F_{4}\frac{\unicode[STIX]{x2202}^{2}u}{\unicode[STIX]{x2202}\unicode[STIX]{x1D709}\unicode[STIX]{x2202}\unicode[STIX]{x1D702}}-F_{5}\frac{\unicode[STIX]{x2202}^{2}u}{\unicode[STIX]{x2202}\unicode[STIX]{x1D709}^{2}}+F_{6}u\frac{\unicode[STIX]{x2202}u}{\unicode[STIX]{x2202}\unicode[STIX]{x1D709}}\\ \displaystyle \hspace{-128.39996pt}+\,F_{6}\frac{\unicode[STIX]{x2202}P}{\unicode[STIX]{x2202}\unicode[STIX]{x1D709}}+F_{2}\frac{\unicode[STIX]{x2202}p}{\unicode[STIX]{x2202}\unicode[STIX]{x1D702}}=0,\\ \displaystyle -\frac{\unicode[STIX]{x2202}^{2}v}{\unicode[STIX]{x2202}\unicode[STIX]{x1D702}^{2}}-F_{1}\frac{\unicode[STIX]{x2202}v}{\unicode[STIX]{x2202}\unicode[STIX]{x1D702}}+F_{2}u\frac{\unicode[STIX]{x2202}v}{\unicode[STIX]{x2202}\unicode[STIX]{x1D702}}+F_{3}v\frac{\unicode[STIX]{x2202}v}{\unicode[STIX]{x2202}\unicode[STIX]{x1D702}}-F_{4}\frac{\unicode[STIX]{x2202}^{2}v}{\unicode[STIX]{x2202}\unicode[STIX]{x1D709}\unicode[STIX]{x2202}\unicode[STIX]{x1D702}}-F_{5}\frac{\unicode[STIX]{x2202}^{2}v}{\unicode[STIX]{x2202}\unicode[STIX]{x1D709}^{2}}+F_{6}u\frac{\unicode[STIX]{x2202}v}{\unicode[STIX]{x2202}\unicode[STIX]{x1D709}}+F_{3}\frac{\unicode[STIX]{x2202}p}{\unicode[STIX]{x2202}\unicode[STIX]{x1D702}}=0\end{array}\right\}\quad & & \displaystyle\end{eqnarray}$$

subject to boundary conditions:

(2.10) $$\begin{eqnarray}\displaystyle \boldsymbol{u}=0\quad \text{at }\unicode[STIX]{x1D702}=\pm 1. & & \displaystyle\end{eqnarray}$$

The coefficients in the above equations are defined as:

(2.11) $$\begin{eqnarray}\displaystyle \left.\begin{array}{@{}c@{}}\displaystyle F_{1}=-\unicode[STIX]{x1D6FC}^{2}\frac{2F_{9}F_{13}-F_{8}F_{10}}{F_{7}},\quad F_{2}=Re\frac{F_{4}}{2\unicode[STIX]{x1D6FC}},\quad F_{3}=2Re\frac{F_{8}}{F_{7}},\\ \displaystyle F_{4}=2\unicode[STIX]{x1D6FC}^{2}\frac{F_{9}F_{8}}{F_{7}},\quad F_{5}=\unicode[STIX]{x1D6FC}^{2}\frac{F_{8}^{2}}{F_{7}},\quad F_{6}=Re\frac{F_{5}}{\unicode[STIX]{x1D6FC}},\\ \displaystyle F_{7}=4+\unicode[STIX]{x1D6FC}^{2}F_{9}^{2},\quad F_{8}=2-2S\cos (\unicode[STIX]{x1D709}),\quad F_{9}=-2S\unicode[STIX]{x1D702}\sin (\unicode[STIX]{x1D709}),\\ \displaystyle F_{10}=-2S\unicode[STIX]{x1D702}\cos (\unicode[STIX]{x1D709}),\quad F_{11}=\unicode[STIX]{x1D6FC}\frac{F_{9}}{F_{8}},\quad F_{12}=\frac{2}{F_{8}},\\ \displaystyle F_{13}=2S\sin (\unicode[STIX]{x1D709}).\end{array}\right\} & & \displaystyle\end{eqnarray}$$

The solution to (2.9) is represented in terms of power expansions of the form:

(2.12) $$\begin{eqnarray}\displaystyle \left.\begin{array}{@{}c@{}}\displaystyle u=u_{0}+\unicode[STIX]{x1D6FC}u_{1}+O(\unicode[STIX]{x1D6FC}^{2}),\\ \displaystyle v=v_{0}+\unicode[STIX]{x1D6FC}v_{1}+O(\unicode[STIX]{x1D6FC}^{2}),\\ \displaystyle p=\unicode[STIX]{x1D6FC}^{-1}p_{-1}+p_{0}+O(\unicode[STIX]{x1D6FC}).\end{array}\right\} & & \displaystyle\end{eqnarray}$$

Substituting (2.12) into (2.9) and retaining the two leading-order terms leads to a sequence of problems with the leading-order system of the form:

(2.13) $$\begin{eqnarray}\displaystyle \left.\begin{array}{@{}c@{}}\displaystyle -\frac{\unicode[STIX]{x2202}^{2}u_{0}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D702}^{2}}+\frac{Re{F_{8}}^{2}}{4}\frac{\unicode[STIX]{x2202}p_{-1}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D709}}=0,\\ \displaystyle \frac{\unicode[STIX]{x2202}p_{-1}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D702}}=0,\\ \displaystyle \frac{\unicode[STIX]{x2202}v_{0}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D702}}=0.\end{array}\right\} & & \displaystyle\end{eqnarray}$$

The next-order system has the form:

(2.14) $$\begin{eqnarray}\displaystyle \left.\begin{array}{@{}c@{}}\displaystyle -\frac{\unicode[STIX]{x2202}^{2}u_{1}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D702}^{2}}+\frac{ReF_{9}F_{8}}{4}u_{0}\frac{\unicode[STIX]{x2202}u_{0}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D702}}+\frac{ReF_{8}}{2}v_{1}\frac{\unicode[STIX]{x2202}u_{0}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D702}}+\frac{ReF_{8}^{2}}{4}u_{0}\frac{\unicode[STIX]{x2202}u_{0}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D709}}+\frac{ReF_{8}^{2}}{4}\frac{p_{0}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D709}}=0,\\ \displaystyle \frac{\unicode[STIX]{x2202}p_{0}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D702}}=0,\\ \displaystyle \frac{\unicode[STIX]{x2202}u_{0}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D709}}+\frac{F_{9}}{F_{8}}\frac{\unicode[STIX]{x2202}u_{0}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D702}}+\frac{2}{F_{8}}\frac{\unicode[STIX]{x2202}v_{1}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D702}}=0.\end{array}\right\} & & \displaystyle\end{eqnarray}$$

The solution of (2.13) is:

(2.15) $$\begin{eqnarray}\displaystyle \left.\begin{array}{@{}c@{}}\displaystyle u_{0}=2M\frac{1-\unicode[STIX]{x1D702}^{2}}{F_{8}},\\ \displaystyle v_{0}=0,\\ \displaystyle p_{-1}=-\frac{H(2+S^{2})}{\sqrt{(1-S^{2})^{5}}}\arctan \left[\frac{\sqrt{1-S}}{\sqrt{1+S}}\tan \left(\frac{\unicode[STIX]{x1D709}}{2}\right)-\frac{\unicode[STIX]{x03C0}}{2}\right]\\ \displaystyle \hspace{31.20007pt}-\,\frac{0.5HS\sin (\unicode[STIX]{x1D709}-\unicode[STIX]{x03C0})[-3S\cos (\unicode[STIX]{x1D709}-\unicode[STIX]{x03C0})+S^{2}-4]}{(1-S^{2})^{2}[1+S\cos (\unicode[STIX]{x1D709}-\unicode[STIX]{x03C0})]^{2}}\\ \displaystyle \hspace{-79.20007pt}-\,\frac{H\unicode[STIX]{x03C0}(2+S^{2})}{\sqrt{(1-S^{2})^{5}}}\left\lfloor \frac{\unicode[STIX]{x1D709}}{2\unicode[STIX]{x03C0}}\right\rfloor \end{array}\right\} & & \displaystyle\end{eqnarray}$$

and the solution of (2.14) is:

(2.16) $$\begin{eqnarray}\displaystyle \left.\begin{array}{@{}c@{}}\displaystyle u_{1}=ReM^{2}\frac{F_{13}}{F_{8}}\left(-\frac{1}{30}\unicode[STIX]{x1D702}^{6}+\frac{1}{6}\unicode[STIX]{x1D702}^{4}-\frac{11}{70}\unicode[STIX]{x1D702}^{2}+\frac{5}{210}\right),\\ \displaystyle v_{1}=\frac{M}{F_{8}}[-S\sin (\unicode[STIX]{x1D709})(\unicode[STIX]{x1D702}^{3}+\unicode[STIX]{x1D702}^{2}-\unicode[STIX]{x1D702}-1)+S\sin (\unicode[STIX]{x1D709})(-\unicode[STIX]{x1D702}^{3}+\unicode[STIX]{x1D702}^{2}+\unicode[STIX]{x1D702}-1)],\\ \displaystyle p_{0}=-\frac{48M^{2}}{35F_{8}^{2}}+c_{0},\end{array}\right\} & & \displaystyle\end{eqnarray}$$

where $M=\sqrt{(1-S^{2})^{5}}/(1+S^{2}/2)$ , $H=2M/Re$ , $c_{0}=12M/(35(1-S)^{2})-(p_{-1}(0))/\unicode[STIX]{x1D6FC}$ is the pressure normalization such that $p=0$ at $\unicode[STIX]{x1D709}=0$ . The last operator used in the expression for $p_{-1}$ is the floor function $\lfloor \cdot \rfloor$ which is added to remove spurious discontinuities associated with the inverse tangent function (Jeffrey & Rich Reference Jeffrey and Rich1994; Mohammadi & Floryan Reference Mohammadi and Floryan2012). The flow rate calculated using the approximate solution (2.15)–(2.16), is:

(2.17) $$\begin{eqnarray}\displaystyle Q=\frac{4}{3}\frac{\sqrt{(1-S^{2})^{5}}}{1+\displaystyle \frac{S^{2}}{2}}+O(\unicode[STIX]{x1D6FC}^{2}), & & \displaystyle\end{eqnarray}$$

i.e. it is not equal to the reference value of $Q_{r}=4/3$ for $\unicode[STIX]{x1D6FC}\rightarrow 0$ , unless $S=0$ . One can define a Reynolds number based on the flow rate in the form of

(2.18) $$\begin{eqnarray}\displaystyle Re_{av}=\frac{Q}{2\unicode[STIX]{x1D708}} & & \displaystyle\end{eqnarray}$$

as it provides a useful assessment tool for the state of the flow. The reference value of $Re_{av}$ for the smooth channel relates to $Re$ defined in § 2.1 as:

(2.19) $$\begin{eqnarray}\displaystyle \frac{Re_{av}}{Re}=\frac{6}{4}. & & \displaystyle\end{eqnarray}$$

2.2.2 Corrugations with wavelength $O(1)$

The analysis of flows in channels with a long-wavelength convergence–divergence discussed in the previous section demonstrates the existence of stationary states and the absence of flow separation. Reduction of the wavelength leads to the onset of separation, as illustrated in figure 2 for $\unicode[STIX]{x1D6FC}=1$ . This separation occurs in the expanding section of the channel, and is stationary and symmetric in the $y$ -direction when $S$ is small enough (see figure 2 a). An increase of $S$ while keeping the same $Re$ results in a transition to an oscillatory separation which does not preserve the symmetry (see figure 2 b). Further reduction of the wavelength leads to an analogous sequence of changes in the flow topology, as illustrated in figure 3 for $\unicode[STIX]{x1D6FC}=3$ . Steady separation occurs for $S=0.075$ (see figure 3 a), the size of the separation zone increases when measured in terms of the fraction of the trough height as $S$ increases (see figure 3 b), it almost completely fills in the trough for $S=0.2$ , pushing the stream away from the walls (see figure 3 c), and it becomes oscillatory for $S=0.225$ (see figure 3 d). At $\unicode[STIX]{x1D6FC}=8$ (not shown) no oscillatory separation has been found in the range of $S$ and $Re$ considered in this study.

Figure 2.

Streamlines of the stationary ( $S=0.125$ , $\unicode[STIX]{x1D6EF}=5.48\times 10^{-6}$ ) (a) and a snapshot of an oscillatory ( $S=0.150$ , $\unicode[STIX]{x1D6EF}=1.15\times 10^{-1}$ ) (b) velocity field. $\unicode[STIX]{x1D6FC}=1$ , $Re=1500$ . Computational boxes containing one corrugation wavelength have been used to produce these results.

Figure 3.

Development of the recirculation zone for $\unicode[STIX]{x1D6FC}=3$ , $Re=2300$ , at different corrugation amplitudes: (a) $S=0.075$ , (b) $S=0.150$ , (c) $S=0.2$ and (d) $S=0.225$ (a snapshot of an oscillatory solution). Computational boxes containing one corrugation wavelength have been used to produce these results.

The above discussion demonstrates the existence of three states, i.e. the stationary states with no separation, the stationary states with symmetric separation and the oscillatory states with non-symmetric separation, and leads to the question of the determination of conditions leading to the appearance of these states. Steady separation can be easily detected through an analysis of the flow field. Identification of the onset of the oscillatory separation is more complex as the numerical solution relies on time simulations and stationary states are detected by examining the time history of the relevant flow properties, e.g. kinetic energy, discussed in § 2.1. The convergence of the numerical process is always oscillatory and, thus, one needs to develop a reliable criterion for separation of the purely numerical oscillations from the physical oscillations, especially at the onset when the physical oscillations have very small amplitudes. Since it has been observed by Cho et al. (Reference Cho, Kim and Shin1998) that stationary flows through symmetric channels exhibit symmetry with respect to the centreline (see figures 2 a, 3 a,b,c; $u$ is symmetric, $v$ is anti-symmetric, $p$ is symmetric) and transition to the oscillatory flow happens through breaking of this symmetry, we use change in the flow symmetry as a criterion for identification of the onset of flow unsteadiness. Accordingly, we define the symmetricity parameter $\unicode[STIX]{x1D6EF}$ as:

(2.20) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6EF}=\frac{1}{N}\mathop{\sum }_{1}^{N}\int _{0}^{y_{u}(x_{n})}(|u(y)-u(-y)|+|v(y)+v(-y)|+|p(y)-p(-y)|)\,\text{d}y & & \displaystyle\end{eqnarray}$$

over a snap shot of the flow field where $N$ is the number of sample points over one corrugation wavelength and $y_{u}(x_{n})$ denotes point at the upper wall at $x_{n}$ . $N=9$ sample points have been used in this work with all points spaced evenly in the $x$ -direction from $x=0.1\unicode[STIX]{x1D706}$ to $x=0.9\unicode[STIX]{x1D706}$ .

Figure 4.

Variations of the size of the separation zone measured at the widest channel opening as a function of $Re$ for stationary flows.

$S=0.025$ , 

$S=0.05$ , 

$S=0.075$ , 

$S=0.1$ , 

$S=0.125$ , 

$S=0.15$ , 

$S=0.175$ , 

$S=0.2$ . Computational boxes containing one corrugation wavelength have been used to produce these results.

We shall now discuss properties of stationary flows. Here we define such flows as flows with $\unicode[STIX]{x1D6EF}<10^{-3}$ . We shall describe a more precise stationarity measure in the next section. We shall limit our interest to $S<0.2$ and $Re<3500$ . The reasons for selecting this range of parameters will become more obvious as we progress through the discussion. Data displayed in figure 4 provide information about variations of the size of the separation zone as a function of the channel geometry and flow conditions. We shall use the height of the separation bubble $H$ at the widest channel opening as a measure of the size of separation. For long enough corrugations (figure 4 a), it is possible to avoid separation for a given $S$ if $Re$ is small enough; $Re$ needs to reach a certain minimum threshold before separation can take place. An increase of $S$ reduces this threshold. Once separation appears, an increase of $Re$ results in an increase of $H$ with $H/S=2$ being the upper limit corresponding to the bubble completely filling in the trough. This increase is very rapid for large $S$ (see $S=0.2$ in figure 4 a for $\unicode[STIX]{x1D6FC}=2$ ) and very gradual for smaller $S$ (see $S=0.025$ , $\unicode[STIX]{x1D6FC}=2$ in figure 4 a). Reduction of the groove wavelength results in a qualitatively similar process with the bubble size increasing more rapidly than for $\unicode[STIX]{x1D6FC}=2$ (compare results in figure 4 b for $\unicode[STIX]{x1D6FC}=5$ with results for $\unicode[STIX]{x1D6FC}=2$ in figure 4 a). Further reduction of the groove wavelength continues this process (see results for $\unicode[STIX]{x1D6FC}=8$ in figure 4 c). It can be seen that in the case of grooves with shorter wavelengths, an increase of $S$ results in the trough being filled very rapidly by the separation bubble (see figure 4 b,c) resulting in the stream being pushed away from the walls.

The presence of grooves affects the flow topology, as discussed above, and results in an increase of the flow resistance. As the pressure gradient is kept constant, increase of the resistance manifests itself through the reduction of the flow rate $Q$ . Figure 5 illustrates the character of these variations as a function of the geometry and the flow conditions. Reduction of $Q$ is rather small and insensitive to variations of $Re$ for the long-wavelength, small-amplitude grooves. These variations become much bigger and very sensitive to $Re$ for larger amplitudes (see results figure 5 a for $\unicode[STIX]{x1D6FC}=0.5$ ). The reduction of the groove wavelength to $\unicode[STIX]{x1D6FC}=3$ (figure 5 b) and $\unicode[STIX]{x1D6FC}=8$ (figure 5 c) results in a significant reduction of $Q$ as the separation bubbles fill in a significant portion of the cross-sectional flow area. At the same time the flow losses become insensitive to variations of $Re$ and $\unicode[STIX]{x1D6FC}$ (compare figure 5 b,c). The latter fact is well illustrated in figure 6 displaying variations of $Q$ as a function of $\unicode[STIX]{x1D6FC}$ . The relative sensitivity of these losses for small $\unicode[STIX]{x1D6FC}$ and large $S$ and their relative insensitivity for large $\unicode[STIX]{x1D6FC}$ and small $S$ are well illustrated.

Figure 5.

Variations of the flow rate $Q$ as a function of $Re$ for stationary flows.

$S=0.025$ , 

$S=0.05$ , 

$S=0.075$ , 

$S=0.1$ , 

$S=0.125$ , 

$S=0.15$ , 

$S=0.175$ , 

$S=0.2$ . Computational boxes containing one corrugation wavelength have been used to produce these results.

Figure 6.

Variation of the flow rate $Q$ as a function of the corrugation wavenumber $\unicode[STIX]{x1D6FC}$ for stationary flows. Thick lines identify limits for $\unicode[STIX]{x1D6FC}\rightarrow 0$ . Computational boxes containing one corrugation wavelength have been used to produce these results.

The reduction of $Q$ leads to the slowing down of the flow and reduction of the effective Reynolds number. This reduction is marginal for small $\unicode[STIX]{x1D6FC}$ (see figure 7 a) but becomes more and more significant as $\unicode[STIX]{x1D6FC}$ increases (see figure 7 b,c). These changes in the effective Reynolds number need to be accounted for in the discussion of stability properties of the flow.

Figure 7.

Variations of the flow rate based Reynolds number $Re_{av}$ as a function of $Re$ for stationary flows.

$S=0.025$ , 

$S=0.05$ , 

$S=0.075$ , 

$S=0.1$ , 

$S=0.125$ , 

$S=0.15$ , 

$S=0.175$ , 

$S=0.2$ . The dashed line corresponds to the smooth channel where $Re_{av}=6/4Re$ . Computational boxes containing one corrugation wavelength have been used to produce these results.

2.3 Transition to oscillatory flow

Figure 8 illustrates variations of $\unicode[STIX]{x1D6EF}$ as a function of $Re$ , for a few illustrative geometries. The investigations were performed using computational boxes containing one corrugation wavelength. It is quite clear that there exists a critical value of $Re$ which separates two distinctive states, one being stationary and the other being oscillatory. Variations of the critical $Re$ as a function of $\unicode[STIX]{x1D6FC}$ and $S$ illustrated in figure 9 suggest that oscillatory states occur for smaller $Re$ and $S$ for grooves with smaller $\unicode[STIX]{x1D6FC}$ . Crossing of lines for $\unicode[STIX]{x1D6FC}=1$ and $\unicode[STIX]{x1D6FC}=2$ suggests otherwise, but this is an artefact of the method used in this analysis related to the size of the computational box, as explained later in this section. The consequences of transition to an oscillatory state are illustrated in figure 10 displaying variations of the flow rate $Q$ as a function of $Re$ . A rapid drop in $Q$ accompanies such a transition. A hysteresis in the change of $Q$ has been observed depending on whether $Re$ is increasing and the system transits from a stationary to the oscillatory state, or $Re$ is decreasing and the system transits from the oscillatory to the stationary state.

Figure 8.

Variations of the symmetricity parameter (see (2.20) for definition) as a function of $Re$ . Flows with $\unicode[STIX]{x1D6EF}>10^{-3}$ are considered to be non-stationary. Internal boxes display instantaneous pressure fields at the jump points. Computational boxes containing one corrugation wavelength have been used to produce these results.

Figure 9.

Borders between the stationary and oscillatory states in the $(Re,S)$ -plane identified using $\unicode[STIX]{x1D6EF}=10^{-3}$ criterion. Locations of the borders have been determined approximately due to the subcritical character of transition between both states. The intersection of curves for $\unicode[STIX]{x1D6FC}=1$ and $\unicode[STIX]{x1D6FC}=2$ results from the change of the wavenumber $\unicode[STIX]{x1D6FF}$ of the most unstable wave as a function of the amplitude $S$ . Computational boxes containing one corrugation wavelength have been used to produce these results.

Figure 10.

Variations of the flow rate $Q$ as a function of $Re$ for $\unicode[STIX]{x1D6FC}=1$ and $S=0.2$ . $\unicode[STIX]{x1D6EF}=10^{-3}$ criterion has been used to distinguish between the stationary and oscillatory states. Computational boxes containing one corrugation wavelength have been used to produce these results.

The features of transition between the stationary and oscillatory states described above suggest that the mechanics of the process can be described using stability techniques and the onset of the oscillatory state represents a subcritical Hopf bifurcation. It has been argued that the instability has the form of a two-dimensional travelling wave analogous to the TS waves (Ghaddar et al. Reference Ghaddar, Korczak, Mikic and Patera1986; Guzmȧn & Amon Reference Guzmȧn and Amon1994; Cho et al. Reference Cho, Kim and Shin1998) arising at Reynolds numbers significantly smaller than for the flow in a smooth channel.

In order to determine the critical conditions leading to transition from stationary to oscillatory states, we represent the total flow quantities as a superposition of the stationary flow, described in § 2.2.2 and a small disturbance, i.e.

(2.21) $$\begin{eqnarray}\displaystyle \left.\begin{array}{@{}c@{}}\displaystyle \boldsymbol{U}_{T}(x,y,z,t)=\boldsymbol{U}(x,y)+\boldsymbol{u}_{p}(x,y,z,t),\\ \displaystyle P_{T}(x,y,z,t)=P(x,y)+p_{p}(x,y,z,t),\end{array}\right\} & & \displaystyle\end{eqnarray}$$

where the subscripts $T$ and $p$ refer to the total flow quantities and small disturbances, respectively. Substitution of (2.21) into the field equation (2.2) and linearisation lead to the disturbance equation of the form

(2.22) $$\begin{eqnarray}\displaystyle \left.\begin{array}{@{}c@{}}\displaystyle \unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}_{\boldsymbol{p}}=0,\\ \displaystyle \frac{\unicode[STIX]{x2202}u_{p}}{\unicode[STIX]{x2202}t}+\boldsymbol{u}_{\boldsymbol{p}}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\boldsymbol{U}+\boldsymbol{U}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\boldsymbol{u}_{\boldsymbol{p}}=-\frac{1}{\unicode[STIX]{x1D71A}}\unicode[STIX]{x1D735}p_{p}+\frac{1}{Re}\unicode[STIX]{x1D6FB}^{2}\boldsymbol{u}_{\boldsymbol{ p}}.\end{array}\right\} & & \displaystyle\end{eqnarray}$$

Disturbances are expressed as

(2.23) $$\begin{eqnarray}\displaystyle q_{p}(x,y,z,t)=Q_{p}(x,y,z)\text{e}^{\text{i}\unicode[STIX]{x1D70E}t}+\text{c.c.}, & & \displaystyle\end{eqnarray}$$

where $q_{p}$ stands for any disturbance quantity, $\unicode[STIX]{x1D70E}=\unicode[STIX]{x1D70E}_{r}+\text{i}\unicode[STIX]{x1D70E}_{i}$ is a complex amplification rate whose real part stands for the frequency and imaginary part represents the amplification rate, $Q_{p}(x,y,z)$ is an $x$ - and $z$ -periodic amplitude function and c.c. stands for the complex conjugate. Substitution of (2.23) into (2.22) leads to an eigenvalue problem for the partial differential equation describing $Q_{p}$ . Discretization of this equation results in an algebraic eigenvalue problem for $\unicode[STIX]{x1D70E}$ of the form:

(2.24) $$\begin{eqnarray}\displaystyle Au_{j}=\unicode[STIX]{x1D70E}_{j}u_{j}, & & \displaystyle\end{eqnarray}$$

which is solved numerically.

It has been shown that the unstable disturbances are equivalent to the TS waves in the limit of a smooth channel (Floryan & Floryan Reference Floryan and Floryan2010). Cho et al. (Reference Cho, Kim and Shin1998) indicates that the three-dimensional oscillatory mode can be interpreted as an oblique (three-dimensional) TS wave. It is convenient for the further discussion to explicitly extract the wave component from the functional disturbance representation resulting in the following disturbance form:

(2.25) $$\begin{eqnarray}\displaystyle q_{p}(x,y,z,t)=\tilde{Q}_{p}(x,y)\text{e}^{\text{i}(\unicode[STIX]{x1D6FF}x+\unicode[STIX]{x1D707}z-\unicode[STIX]{x1D70E}t)}+\text{c.c.}, & & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D6FF}$ and $\unicode[STIX]{x1D707}$ are the real wavenumbers in the $x$ - and $z$ -directions, respectively and $\tilde{Q}_{p}(x,y)$ stands for the $x$ -periodic amplitude function. As it has been reported by Cho et al. (Reference Cho, Kim and Shin1998) and Floryan & Floryan (Reference Floryan and Floryan2010) that the two-dimensional waves are less stable (having lower critical Reynolds number) than the three-dimensional waves, the further discussion is limited to such waves and, accordingly, we set $\unicode[STIX]{x1D707}=0$ . We shall return to this issue later in the discussion.

Discretization leading to (2.24) is based on the same spectral finite-element method as used in § 2.1 and uses periodic boundary conditions in the $x$ -direction. As a result, the $x$ -size of the computational box determines admissible $\unicode[STIX]{x1D6FF}$ values (Blancher et al. Reference Blancher, Le Guer and El Omari2015).

Figure 11.

Variations of the critical Reynolds numbers $Re_{c}$ of the two most unstable, two-dimensional travelling waves (a) and the corresponding critical frequencies $\unicode[STIX]{x1D70E}_{r,c}$ (b) as functions of the groove amplitude $S$ for $\unicode[STIX]{x1D6FC}=1$ . Computational boxes containing one corrugation wavelength have been used to produce these results. $\unicode[STIX]{x1D6FF}=\unicode[STIX]{x1D6FF}_{cr}$ indicates the critical Reynolds number calculated using the immersed boundary conditions (IBC) method (Floryan & Floryan Reference Floryan and Floryan2010) which permits a nearly continuous variations of $\unicode[STIX]{x1D6FF}$ .

Assuming that the length of the computational box is $L=n\ast \unicode[STIX]{x1D706}$ , where $n$ is the number of corrugations contained in the box, the discretized operator can only capture waves with wavelengths $L/m$ where $m$ is an integer. This limits the disturbance wavelengths which can be captured using this particular box to $\unicode[STIX]{x1D6FF}=\unicode[STIX]{x1D6FC}m/n$ , i.e. it is not possible to vary $\unicode[STIX]{x1D6FF}$ values in a continuous manner. Nevertheless, useful conclusions regarding characteristics of various waves can be extracted by properly adjusting the length of the computational box. In the analysis, boxes with length $L=2\unicode[STIX]{x03C0}$ were used to study waves with $\unicode[STIX]{x1D6FF}=1,2,\ldots ,m$ . In the case of $\unicode[STIX]{x1D6FC}=3$ , boxes containing $n=2$ and $5$ corrugations were used to study waves with non-integer $\unicode[STIX]{x1D6FF}$ values. Verification of the accuracy of this method is discussed in appendix B. It is known that flow in a smooth channel becomes unstable at $Re_{cr}=5772$ and the two-dimensional (TS) wave with the wavenumber $\unicode[STIX]{x1D6FF}_{c}=1.02$ travelling downstream plays the critical role. Floryan & Floryan (Reference Floryan and Floryan2010) reports that in the case of small groove amplitudes the critical wavenumber changes marginally. Figure 11(a) illustrates the dependence of the critical Reynolds number on the amplitude $S$ determined for $\unicode[STIX]{x1D6FC}=1$ using a computational box with length equal to one groove wavelength. Variations of frequencies $\unicode[STIX]{x1D70E}_{r}$ of the two dominant instability modes are displayed in figure 11(b). It can be seen that the wave with frequency $\unicode[STIX]{x1D70E}_{r}\approx 0.2$ plays the critical role until $S$ reaches a certain critical magnitude and then the wave with frequency $\unicode[STIX]{x1D70E}_{r}\approx 0.6$ begins to play the critical role. Tracking of frequencies of these two waves in the limit $S\rightarrow 0$ demonstrates that they connect to TS waves with the wavenumbers $\unicode[STIX]{x1D6FF}=1$ and $\unicode[STIX]{x1D6FF}=2$ (figure 12). These results suggest that growth of the groove amplitude results in the reduction of the effective channel opening giving preference to waves with shorter wavelengths. Figure 13 illustrates the behaviour of these waves in channels with different corrugation wavenumbers. The critical Reynolds numbers for both waves in channels with long wavelength grooves ( $\unicode[STIX]{x1D6FC}=0.5$ ; see figure 13 a) decrease monotonically as $S$ increases. The wave with $\unicode[STIX]{x1D6FF}=2$ becomes the dominant one at $S\succsim 0.07$ . When the groove wavelength is reduced to $\unicode[STIX]{x1D6FC}=5$ (see figure 13 b), $Re_{cr}$ for wave with $\unicode[STIX]{x1D6FF}=1$ reaches a minimum for $S\approx 0.1$ and increases for larger $S$ . The wave with $\unicode[STIX]{x1D6FF}=2$ begins to dominate for $S\succsim 0.12$ . When the groove wavelength is further reduced (see figure 13 for $\unicode[STIX]{x1D6FC}=8$ ), the non-monotonic variations for $Re_{cr}$ for the wave with $\unicode[STIX]{x1D6FF}=1$ are more pronounced and amplitude $S$ has to reach a higher value before the wave with $\unicode[STIX]{x1D6FF}=2$ becomes critical. In this particular case, i.e.  $\unicode[STIX]{x1D6FC}=8$ , the wave with $\unicode[STIX]{x1D6FF}=2$ dominates for $S\succsim 0.17$ .

Figure 12.

Variations of the amplification rate $\unicode[STIX]{x1D70E}_{i}$ (a), and the frequency $\unicode[STIX]{x1D70E}_{r}$ (b) of the travelling wave instability for $\unicode[STIX]{x1D6FC}=1$ for the two most unstable two-dimensional waves as functions of the groove amplitude $S$ for $Re=1000$ . Computational boxes containing one corrugation wavelength have been used to produce these results.

Figure 13.

Variations of the critical conditions for the two most unstable, two-dimensional waves as functions of the groove amplitude $S$ : (a–c) display the critical Reynolds numbers $Re_{c}$ for the groove wavenumbers $\unicode[STIX]{x1D6FC}=0.5,5,8$ , respectively, and (d–f) display the corresponding critical frequencies $\unicode[STIX]{x1D70E}_{r,c}$ . Computational boxes of length $2\unicode[STIX]{x03C0}$ have been used to produce these results.

The above results were produced using computational boxes with lengths selected to contain integer numbers of corrugation and waves‘ wavelengths. One should expect continuous variations of the critical $\unicode[STIX]{x1D6FF}$ as a function of $S$ but this process cannot be followed using the present numerical method. The continuous variations of $\unicode[STIX]{x1D6FF}$ are illustrated in qualitative sense in figure 14 for the groove wavenumber $\unicode[STIX]{x1D6FC}=3$ . The lengths of the computational boxes have been adjusted to allow for capturing waves with different $\unicode[STIX]{x1D6FF}$ values. It can be seen that an increase of $S$ gives preference to waves with shorter wavelengths.

Figure 14.

Variations of the critical Reynolds number $Re_{c}$ (a) and the critical frequency $\unicode[STIX]{x1D70E}_{r,c}$ (b) for the least stable two-dimensional travelling waves for the corrugation wavenumber $\unicode[STIX]{x1D6FC}=3$ as functions of the corrugation amplitude $S$ . Computational boxes with the properly selected lengths have been used to study waves with fractional $\unicode[STIX]{x1D6FF}$ values.

3.1 Three-dimensional travelling wave

It is known that the two-dimensional waves are less attenuated than the three-dimensional waves in channels with small amplitude of convergence–divergence $S$ (Cho et al. Reference Cho, Kim and Shin1998; Floryan & Floryan Reference Floryan and Floryan2010). It remains to be shown that this relation holds when the amplitude $S$ is increased. The solution procedure is identical to that described in § 2.3 with small changes explained below. The disturbances are assumed to be in the form:

(3.1) $$\begin{eqnarray}\displaystyle q_{p}(x,y,z,t)=\tilde{\dot{Q}}_{p}(x,y)\text{e}^{\text{i}(\unicode[STIX]{x1D6FF}x+\unicode[STIX]{x1D707}z-\unicode[STIX]{x1D70E}t)}+\text{c.c.} & & \displaystyle\end{eqnarray}$$

resulting in an eigenvalue problem for the partial differential equations describing $\tilde{\dot{Q}}_{p}(x,y)$ . One only needs to discretize $Q_{p}$ in the $(x,y)$ -plane. The discretization process is the same as described in § 2.1 resulting in an algebraic eigenvalue problem where $\unicode[STIX]{x1D707}$ and $\unicode[STIX]{x1D6FF}$ play the role of parameters. Results presented in table 1 demonstrate an increase of the critical Reynolds number as the wave inclination angle increases for small and large amplitudes of corrugation $S$ . It can be concluded that it is unlikely that the flow three-dimensionalization can be caused by the three-dimensional travelling waves.

Table 1.

Variations of the critical Reynolds number $Re_{cr}$ for oblique travelling waves with the wave vector ( $\unicode[STIX]{x1D6FF}=1,2$ , $\unicode[STIX]{x1D707}$ ) as a function of the wave inclination angle expressed in terms of the spanwise wavenumber $\unicode[STIX]{x1D707}$ for grooves with $\unicode[STIX]{x1D6FC}=1$ .

3.2 Streamwise vortices

We shall now focus our attention on the streamwise vortices. Channel modulations force a change of direction of the flow creating a centrifugal-force field which, under certain conditions, may lead to the formation of a secondary flow in the form of streamwise vortices. As the centrifugal force drives the instability and viscous effects oppose it, the vortices may appear only when the former one dominates, i.e. when the streamline curvature is large enough to generate centrifugal forces strong enough to overcome friction. It is rather obvious that an increase of $Re$ will lead to the flow destabilization as larger $Re$ signals the reduction of viscous forces. The effects of geometry changes are more difficult to predict without detailed analysis. An initial increase of the groove amplitude increases streamline curvature but an excessive increase leads to flow separation which decreases the curvature. It goes without saying that streamlines of the main stream, and not streamlines inside the separation bubbles, are of consequence.

The conditions leading to the onset of streamwise vortices can be predicted using stability theory. Disturbance quantities are represented using equation (2.23) with $\unicode[STIX]{x1D70E}_{r}=0$ since vortices do not propagate, as opposed to the waves. Rather than solving an eigenvalue problem for the disturbance equations, we solve the complete time-dependent problem described by (2.2), where the flow quantities are represented by a Fourier expansion truncated to $N$ leading modes, i.e.

(3.2) $$\begin{eqnarray}\displaystyle \boldsymbol{u}(x,y,z,t)=\mathop{\sum }_{k=-N}^{k=N}\hat{\boldsymbol{u}}_{k}(x,y,t)\text{e}^{\text{i}k\unicode[STIX]{x1D707}z}, & & \displaystyle\end{eqnarray}$$

where $\hat{\boldsymbol{u}}_{k}=\hat{\boldsymbol{u}}_{-k}^{\ast }$ are the reality conditions and stars identify the complex conjugate. The initial conditions consist of a superposition of the two-dimensional solution, determined in § 2.2.2, and a suitable three-dimensional perturbation. The evolution of the perturbations is followed in time by computing the kinetic energy of the Fourier mode of interest defined as

(3.3) $$\begin{eqnarray}\displaystyle E_{k}(t)=\frac{1}{2}\int _{\unicode[STIX]{x1D6FA}}\Vert \hat{\boldsymbol{u}}_{k}(t)\Vert ^{2}\,\text{d}\unicode[STIX]{x1D6FA}, & & \displaystyle\end{eqnarray}$$

where $k$ denots the mode number. According to (2.23), the energy is expected to vary exponentially, i.e.

(3.4) $$\begin{eqnarray}\displaystyle E_{k}(t)=A_{k}\cdot \text{e}^{B_{k}\cdot t}, & & \displaystyle\end{eqnarray}$$

where $A_{k}$ and $B_{k}$ are constants to be determined. The reader may note that the growth rate for mode $k$ can be expressed as $\unicode[STIX]{x1D70E}_{i,k}=B_{k}/2$ . The modal functions are discretized on a box of length $\unicode[STIX]{x1D706}=2\unicode[STIX]{x03C0}/\unicode[STIX]{x1D6FC}$ using the same spectral finite-element discretization method as used in § 2.2.2. The same temporal discretization as in § 2.2.2 is used to determine the time evolution.

All modes can be traced during time simulation. For each mode a nonlinear least-square fitting of (3.4) to a sample of data is used to determine constants $A_{k}$ and $B_{k}$ as illustrated in figure 15. Since the first mode is the most unstable (Floryan Reference Floryan2003), and to limit the computational cost, only the zeroth mode (mean flow representation) and the first mode (the dominant disturbance mode) have been traced during the majority of simulations. This limits the results to the linear stability theory. We shall discuss higher modes when discussing saturation states.

Figure 15.

Evolution of the kinetic energy of the first Fourier mode (see (3.3)) for flows with $Re=500$ and $Re=2500$ in a channel with $\unicode[STIX]{x1D6FC}=3$ , $S=0.05$ , $\unicode[STIX]{x1D707}=3.5$ . Solid lines illustrate the complete time history while the thick lines identify the data segments used for determination of $A_{k}$ and $B_{k}$ (see (3.4)). Computational box containing one corrugation wavelength has been used to produce these results.

Figure 15 illustrates the data fitting process used to determine $A_{k}$ and $B_{k}$ . Variations of the kinetic energy are plotted as a function of time and the curve fitting produces approximations for the growth rate. Accuracy verification of this procedure is discussed in appendix C.

Figure 16 illustrates variations of the amplification rate as a function of the vortex wavenumber $\unicode[STIX]{x1D707}$ for corrugations with the amplitude $S=0.1$ and with different wavelengths for $Re=1300$ . It can be seen that vortices which are either too large or too small are being attenuated. The range of unstable $\unicode[STIX]{x1D707}$ is a strong function of $\unicode[STIX]{x1D6FC}$ with grooves which are either too long or too short unable to cause the instability. In the former case, the streamline curvature is dictated by the slope of the bounding walls (see figure 2 a). Since this slope is very small, the centrifugal force is too weak to cause the instability. In the latter case, the stream lifts above the grooves (see figure 3) resulting in nearly rectilinear streamlines and a centrifugal force which is also too weak to create the secondary flow.

Figure 16.

Variations of the vortex amplification rate as a function of the vortex wavenumber $\unicode[STIX]{x1D707}$ for different corrugation wavelengths for $Re=1300$ and $S=0.1$ . Computational boxes containing one corrugation wavelength have been used to produce these results.

Figure 17 illustrates the neutral curves in the $(Re,\unicode[STIX]{x1D707})$ -plane for corrugations with the fixed wavelength $\unicode[STIX]{x1D706}=2\unicode[STIX]{x03C0}/\unicode[STIX]{x1D6FC}$ but with the variable amplitude $S$ . The minima of these curves determine the critical conditions, i.e. the critical Reynolds numbers $Re_{cr}$ and the critical vortex wavenumbers $\unicode[STIX]{x1D707}_{cr}$ . It can be seen that an increase of $S$ results in negligible variations of $\unicode[STIX]{x1D707}_{cr}$ but changes of $Re_{cr}$ are more complex and depend on $\unicode[STIX]{x1D6FC}$ . In the case of long-wavelength grooves ( $\unicode[STIX]{x1D6FC}=1$ , figure 17 a), the initial increase of $S$ results in a rapid decrease of $Re_{cr}$ which eventually saturates with the neutral curves for $S=0.15$ and $S=0.20$ nearly overlapping. In the case of the grooves with shorter wavelength ( $\unicode[STIX]{x1D6FC}=3$ , figure 17 b), the initial rapid decrease of $Re_{cr}$ is followed by an increase once $S$ becomes large enough. In this case, most of the neutral curve for $S=0.175$ is above the neutral curve for $S=0.15$ . Figure 17(a) also displays the neutral curves for the travelling wave instability. It can be seen that the two-dimensional waves dictate the critical conditions with $Re_{cr}$ for the waves being larger than $Re_{cr}$ for the vortices.

Figure 17.

The neutral stability curves in the $(Re,\unicode[STIX]{x1D707})$ -plane for different groove amplitudes $S$ . Solid lines correspond to the critical Reynolds numbers for the vortex instability. Dashed lines correspond to the critical Reynolds numbers for the travelling wave instability. ○ Identifies the two-dimensional wave for $S=0.02$ and ▫ identifies such wave for $S=0.2$ . Computational boxes containing one corrugation wavelength have been used to produce these results.

Effects of the groove amplitude are explicitly illustrated in figure 18, which displays variations of the amplification rate for vortices with $\unicode[STIX]{x1D707}_{cr}$ as a function of $S$ for a fixed $Re$ . In all cases the groove amplitude needs to reach a certain minimum before a vortex instability may be initiated. If the groove amplitude is too small, the system response is dominated by the travelling wave instability (Floryan & Floryan Reference Floryan and Floryan2010). In the case of the long-wavelength grooves ( $\unicode[STIX]{x1D6FC}=1$ , figure 18 a) the vortex amplification rate monotonically increases with $S$ up to the maximum $S$ considered in this analysis. The vortices can be amplified for $Re$ as low as $Re\approx 400$ . Reduction of the groove wavelength ( $\unicode[STIX]{x1D6FC}=3$ , figure 18 b) changes the system response qualitatively as an increase of $S$ initially leads to an increase of the amplification followed by its decrease and the eventual vortex stabilization even for the highest $Re$ considered. In this case one needs to use at least $Re\approx 700$ in order to create this instability. Further reduction of the groove wavelength ( $\unicode[STIX]{x1D6FC}=5$ , figure 18 c) follows the same trend. One needs to increase the Reynolds number up to $Re\approx 1000$ in order to create a secondary flow and the vortex stabilization takes place when $S$ increases beyond $S\approx 0.12$ . In the case of the short-wavelength grooves ( $\unicode[STIX]{x1D6FC}=8$ , figure 18 d), one needs to use $Re>2000$ to cause the instability, the range of $S$ that can lead to instability is very narrow and the amplification rates are very small. Further increase of $\unicode[STIX]{x1D6FC}$ eliminates the vortex instability altogether in the range of $S$ and $Re$ considered.

Figure 18.

Variations of the vortex amplification rate at $\unicode[STIX]{x1D707}_{cr}$ as a function of the groove amplitude and the flow Reynolds number. Computational boxes containing one corrugation wavelength have been used to produce these results.

Figure 19 displays variations of the critical Reynolds numbers as a function of the corrugation amplitude. A monotonic decrease of $Re_{cr}$ is observed for the longest groove wavelengths tested ( $\unicode[STIX]{x1D6FC}=0.5,1.0$ ). Starting with $\unicode[STIX]{x1D6FC}=2$ , one can identify an amplitude $S$ which corresponds to the minimum critical Reynolds number. We shall refer to this amplitude as the critical amplitude $S_{cr}$ . As the corrugation wavelength further decreases (see curves for $\unicode[STIX]{x1D6FC}=3,5,8$ in figure 19), the critical amplitude decreases to about $S_{cr}\approx 0.06$ and the minimum $Re_{cr}$ increases to approximately $Re_{cr}\approx 1800$ . It is worth noticing that further decreases of groove wavelength to $\unicode[STIX]{x1D6FC}=10$ do not affect $S_{cr}$ but increases the minimum $Re_{cr}$ .

Figure 19.

Variations of the critical Reynolds number $Re_{cr}$ as a function of the groove amplitude $S$ . Symbols identify the computed data points. Dotted lines result from the data extrapolation. Computational boxes containing one corrugation wavelength have been used to produce these results.

Figure 20.

Variations of the minimum critical Reynolds number $Re_{cr,min}$ , the corresponding critical corrugation amplitude $S_{cr}$ and the corresponding reduction of the flow rate as function of the corrugation wavenumber $\unicode[STIX]{x1D6FC}$ . Computational boxes containing one corrugation wavelength have been used to produce these results.

The above results provide a basis for the identification of the most effective geometric configuration for the generation of the vortices. Here we view the configuration which can amplify vortices at the smallest $Re$ as being the most effective. Figure 20 displays variations of the minimum critical Reynolds number and the corresponding critical corrugation amplitude $S_{cr}$ , as functions of the corrugation wavenumber. Depending on $\unicode[STIX]{x1D6FC}$ , the minimum $Re_{cr}$ changes from $300$ to $2200$ and $S_{cr}$ changes from $0.22$ to $0.06$ . The cost of creating vortices can be expressed in terms of the flow rate reduction caused by the grooves. The flow rate reduction corresponding to the use of the optimal configuration varies as a function of $\unicode[STIX]{x1D6FC}$ from $30\,\%$ to $11\,\%$ (see figure 20).

Figure 21.

Variations of the critical Reynolds number $Re_{cr}$ for the onset of the vortex instability (marked with

) and the two-dimensional travelling wave instability (marked with

) with various wavenumbers $\unicode[STIX]{x1D6FF}$ as a function of the groove amplitude $S$ . Thick solid line without any markers in the insert in (b) has been obtained using IBC method. Data for the vortex instability have been determined using computational boxes containing a single corrugation wavelength. Data for the travelling wave instability have been determined using computational boxes containing different number of corrugation wavelengths; boxes with one corrugation have been used for $\unicode[STIX]{x1D6FC}=0.5$ and $1$ , boxes with five corrugation wavelengths have been used for $\unicode[STIX]{x1D6FC}=5$ , and boxes with two, three and five wavelengths have been used for $\unicode[STIX]{x1D6FC}=3$ .

3.3 Competition between the vortices and the travelling waves

The formation of the streamwise vortices is of the main interest in many applications as it results in better flow mixing through creation of transverse fluid motion. The process of vortex creation can be interfered with by the travelling waves. It is thus important to establish conditions where such vortices are likely to dominate the system dynamics. Figure 21 displays variations of $Re_{cr}$ as a function of the groove amplitude required for the onset of the travelling waves and for the onset of the vortices. In the case of the long-wavelength grooves ( $\unicode[STIX]{x1D6FC}=0.5$ , figure 21 a) the vortices dominate in the whole range of corrugation amplitudes being considered in this study with $Re_{cr}$ for both instabilities monotonically decreasing as $S$ increases, and with $Re_{cr}$ for the vortices being always smaller than $Re_{cr}$ for the waves. The reader may recall from the discussion in § 2.3 that our method of analysis does not allow for continuous variations of the wavenumber for the travelling wave instability and, thus, the results presented in figure 21 for such waves should be viewed as representing spot checks. The differences between the corresponding $Re_{cr}$ and $Re_{cr}$ for the vortices are large enough that one can nevertheless conclude that the waves are irrelevant in channels with such corrugations. Reduction of the groove wavelength to $\unicode[STIX]{x1D6FC}=1$ (see figure 21 b) leads to a wider separation of $Re_{cr}$ with the vortices being more dominant. Further decrease of the groove wavelength to $\unicode[STIX]{x1D6FC}=3$ (see figure 21 c) affects very little variations of $Re_{cr}$ for the waves but changes the character of variations of $Re_{cr}$ for the vortices. $Re_{cr}$ in the latter cases initially decreases and then it starts to increase. As a result, the vortices dominate until $S\approx 0.17$ and the waves take over for larger $S$ . It is unlikely that the vortices can be formed in channels with such large $S$ . In the case of the short-wavelength grooves ( $\unicode[STIX]{x1D6FC}=5$ , figure 21 d) formation of the waves requires the use of $Re$ much larger than in the case of $\unicode[STIX]{x1D6FC}=3$ while the vortices can be created with $Re$ only somewhat larger but only as long as the groove amplitude is not too large. Increase of $S$ above 0.12 results in a rapid stabilization of the vortices and the dominance of the waves. The range of $S$ which permits the formation of the vortices without interference from the waves is much smaller in such conduits.