2.1 Navier–Stokes model

The incompressible Navier–Stokes equations in a smooth plane domain $\unicode[STIX]{x1D6FA}$ read

(2.1a ) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x2202}_{t}\boldsymbol{u}+(\boldsymbol{u}\boldsymbol{\cdot }\unicode[STIX]{x1D735})\boldsymbol{u}=\unicode[STIX]{x1D708}\unicode[STIX]{x0394}\boldsymbol{u}-\unicode[STIX]{x1D735}p, & \displaystyle\end{eqnarray}$$

(2.1b ) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}=0, & \displaystyle\end{eqnarray}$$

where $\boldsymbol{u}(\boldsymbol{x},t)$ is the velocity field, $p$ is the pressure field, $\unicode[STIX]{x1D708}$ is the kinematic viscosity, and we shall denote by $(u,v)$ the two components of $\boldsymbol{u}$ . In order to make formulas a little more concise, we shall in the following often omit to write the time variable explicitly, except when doing so would create an ambiguity.

As spatial domain, we choose

(2.1c ) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6FA}=\{(x,y)\mid x\in \mathbb{T},~y\in ]0,1[\,\}, & & \displaystyle\end{eqnarray}$$

where $\mathbb{T}=\mathbb{R}/\mathbb{Z}$ is the unit circle, which models a periodic channel. Dirichlet boundary conditions are imposed at $y=0$ and $y=1$ ,

(2.1d ) $$\begin{eqnarray}\displaystyle \boldsymbol{u}(x,0,t)=\boldsymbol{u}(x,1,t)=0, & & \displaystyle\end{eqnarray}$$

and we specify an initial condition

(2.1e ) $$\begin{eqnarray}\displaystyle \boldsymbol{u}(x,y,0)=\boldsymbol{u}_{i}(x,y), & & \displaystyle\end{eqnarray}$$

which we shall assume to have zero spatial average. By introducing characteristic velocity and length scales $U$ and $L$ , a Reynolds number can be defined as follows:

(2.2) $$\begin{eqnarray}\displaystyle Re=\frac{UL}{\unicode[STIX]{x1D708}}. & & \displaystyle\end{eqnarray}$$

When discretizing this system, difficulties arise due to the interplay between the divergence condition (2.1b ) and the no-slip boundary condition (2.1d ). We have chosen to work with the vorticity formulation of the NSEs, which eliminates the divergence constraint at the cost of transforming the Dirichlet boundary conditions on $u$ into a non-local integral constraint on $\unicode[STIX]{x1D714}$ . Although they used to be controversial, such formulations are now well established (see Gresho Reference Gresho1991; Weinan & Liu Reference E and Liu1996; Maekawa Reference Maekawa2013), under the condition that the discretization of the integral constraint is properly carried out. Fortunately, our periodic channel geometry allows for an explicit and easy-to-understand approach, which we now present.

The vorticity field $\unicode[STIX]{x1D714}=\unicode[STIX]{x2202}_{x}v-\unicode[STIX]{x2202}_{y}u$ satisfies the transport equation

(2.3) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x2202}_{t}\unicode[STIX]{x1D714}+\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\boldsymbol{u}\unicode[STIX]{x1D714})=\unicode[STIX]{x1D708}\unicode[STIX]{x0394}\unicode[STIX]{x1D714}, & & \displaystyle\end{eqnarray}$$

with initial data

(2.4) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D714}_{i}(x,y):=\unicode[STIX]{x1D735}\times \boldsymbol{u}_{i}(x,y), & & \displaystyle\end{eqnarray}$$

where $\boldsymbol{u}$ is expressed as a function of $\unicode[STIX]{x1D714}$ by means of the streamfunction $\unicode[STIX]{x1D713}$ defined by

(2.5a ) $$\begin{eqnarray}\displaystyle & \displaystyle u=-\unicode[STIX]{x2202}_{y}\unicode[STIX]{x1D713}, & \displaystyle\end{eqnarray}$$

(2.5b ) $$\begin{eqnarray}\displaystyle & \displaystyle v=\unicode[STIX]{x2202}_{x}\unicode[STIX]{x1D713}, & \displaystyle\end{eqnarray}$$

(2.6) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x0394}\unicode[STIX]{x1D713}=\unicode[STIX]{x1D714}. & & \displaystyle\end{eqnarray}$$

From the wall-normal component of (2.1d ) and the fact that $\int _{\unicode[STIX]{x1D6FA}}\boldsymbol{u}=0$ is a constant of motion, a Dirichlet boundary condition for $\unicode[STIX]{x1D713}$ follows,

(2.7) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D713}(x,0,t)=\unicode[STIX]{x1D713}(x,1,t)=0, & & \displaystyle\end{eqnarray}$$

which uniquely determines $\unicode[STIX]{x1D713}$ , and therefore $\boldsymbol{u}$ , as a function of $\unicode[STIX]{x1D714}$ .

To close the problem, the tangential component of (2.1d ), which has not yet been used, needs to be reformulated into the missing boundary condition on $\unicode[STIX]{x1D714}$ necessary because of the presence of a Laplacian in (2.3). A general discussion of this issue has been carried out by Gresho (Reference Gresho1991). In our case, due to the simple geometry, (2.5)–(2.7) can be solved explicitly to get an expression for $\unicode[STIX]{x1D713}$ . For this, we first introduce the Fourier coefficients

(2.8) $$\begin{eqnarray}\displaystyle \widehat{\unicode[STIX]{x1D714}}_{k}(y)=\int _{0}^{1}\unicode[STIX]{x1D714}(x,y)\text{e}^{-\text{i}\unicode[STIX]{x1D705}x}\,\text{d}x, & & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D705}=2\unicode[STIX]{x03C0}k$ , and the corresponding reconstruction formula

(2.9) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D714}(x,y)=\mathop{\sum }_{k\in \mathbb{Z}}\widehat{\unicode[STIX]{x1D714}}_{k}(y)\text{e}^{\text{i}\unicode[STIX]{x1D705}x}, & & \displaystyle\end{eqnarray}$$

which applies similarly for other fields. By (2.6) we then have

(2.10) $$\begin{eqnarray}\displaystyle -\unicode[STIX]{x1D705}^{2}\widehat{\unicode[STIX]{x1D713}}_{k}+\unicode[STIX]{x2202}_{y}^{2}\widehat{\unicode[STIX]{x1D713}}_{k}=\widehat{\unicode[STIX]{x1D714}}_{k}. & & \displaystyle\end{eqnarray}$$

Combining this with the boundary conditions (2.7), we obtain

(2.11a ) $$\begin{eqnarray}\displaystyle \widehat{\unicode[STIX]{x1D713}}_{k}(y) & = & \displaystyle -\frac{1}{2|\unicode[STIX]{x1D705}|(1-\text{e}^{-2|\unicode[STIX]{x1D705}|})}\left((1-\text{e}^{-2|\unicode[STIX]{x1D705}|y})\int _{y}^{1}\widehat{\unicode[STIX]{x1D714}}_{k}(y^{\prime })(\text{e}^{-|\unicode[STIX]{x1D705}|(y^{\prime }-y)}-\text{e}^{-|\unicode[STIX]{x1D705}|(2-y^{\prime }-y)})\,\text{d}y^{\prime }\right.\nonumber\\ \displaystyle & & \displaystyle +\left.(1-\text{e}^{-2|\unicode[STIX]{x1D705}|(1-y)})\int _{0}^{y}\widehat{\unicode[STIX]{x1D714}}_{k}(y^{\prime })(\text{e}^{-|\unicode[STIX]{x1D705}|(y-y^{\prime })}-\text{e}^{-|\unicode[STIX]{x1D705}|(y+y^{\prime })})\,\text{d}y^{\prime }\right),\end{eqnarray}$$

for $k\neq 0$ , and

(2.11b ) $$\begin{eqnarray}\displaystyle \widehat{\unicode[STIX]{x1D713}}_{0}(y)=y\int _{0}^{1}(1-y^{\prime })\widehat{\unicode[STIX]{x1D714}}_{0}(y^{\prime })\,\text{d}y^{\prime }+\int _{0}^{y}(y-y^{\prime })\widehat{\unicode[STIX]{x1D714}}_{0}(y^{\prime })\,\text{d}y^{\prime }. & & \displaystyle\end{eqnarray}$$

Using these expressions, the no-slip boundary condition (2.1d ) can now be reformulated as two linear constraints on $\widehat{\unicode[STIX]{x1D714}}_{k}$ :

(2.12) $$\begin{eqnarray}\displaystyle \forall k,\quad B_{k}^{+}(\widehat{\unicode[STIX]{x1D714}}_{k})=B_{k}^{-}(\widehat{\unicode[STIX]{x1D714}}_{k})=0, & & \displaystyle\end{eqnarray}$$

where

(2.13a,b ) $$\begin{eqnarray}\displaystyle B_{k}^{+}(f)=\int _{0}^{1}\text{e}^{-|\unicode[STIX]{x1D705}|y}f(y)\,\text{d}y,\quad B_{k}^{-}(f)=\int _{0}^{1}\text{e}^{-|\unicode[STIX]{x1D705}|(1-y)}f(y)\,\text{d}y=0, & & \displaystyle\end{eqnarray}$$

for $k\neq 0$ , and

(2.13c,d ) $$\begin{eqnarray}\displaystyle B_{0}^{+}(f)=\int _{0}^{1}yf(y)\,\text{d}y=0\quad \text{and}\quad B_{0}^{-}(f)=\int _{0}^{1}(1-y)f(y)\,\text{d}y=0. & & \displaystyle\end{eqnarray}$$

For numerical purposes, it is better to reformulate these stiff conditions by taking advantage of the diffusion operator, i.e. by applying $B^{+}$ and $B^{-}$ to (2.3). Assuming smooth solutions, the linear constraints commute with the partial time derivative, i.e. we also have $B_{k}^{\pm }(\unicode[STIX]{x2202}_{t}\widehat{\unicode[STIX]{x1D714}}_{k})=0$ . Applying the Fourier transform, for which a similar argument holds, to (2.3) and then applying $B_{k}^{\pm }$ to the resulting equations eliminates the time derivative terms for each wavenumber $k$ and yields

(2.14) $$\begin{eqnarray}\displaystyle B_{k}^{\pm }(\unicode[STIX]{x2202}_{y}^{2}\widehat{\unicode[STIX]{x1D714}}_{k})=\unicode[STIX]{x1D708}^{-1}B_{k}^{\pm }(\widehat{\unicode[STIX]{x1D735}(\boldsymbol{u}\unicode[STIX]{x1D714})}_{k}). & & \displaystyle\end{eqnarray}$$

The above analysis ensures that the system of equations (2.3)–(2.5), (2.11) and either one of (2.12) or (2.14) is equivalent to the original Navier–Stokes system for smooth strong solutions.

2.2 Prandtl–Euler model

We now describe the alternative model for the flow derived by Prandtl (Reference Prandtl1904). Although Prandtl and most later authors used the velocity variable to write down the equations, we present here the equivalent vorticity formulation, since we have found that it leads to a simpler understanding of the phenomena we are interested in.

The starting point is the following ansatz for the vorticity field as $Re\rightarrow \infty$ :

(2.15) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D714}(x,y)=\unicode[STIX]{x1D714}_{E}(x,y)+\unicode[STIX]{x1D708}^{-1/2}\unicode[STIX]{x1D714}_{P}(x,\unicode[STIX]{x1D708}^{-1/2}y)-\unicode[STIX]{x1D708}^{-1/2}\unicode[STIX]{x1D714}_{P}(x,\unicode[STIX]{x1D708}^{-1/2}(1-y))+\unicode[STIX]{x1D714}_{R}(x,y),\quad & & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D714}_{E}(x,y)$ is a smooth function on $\unicode[STIX]{x1D6FA}\times ]0,T[$ , and $\unicode[STIX]{x1D714}_{P}(x,y_{P})$ is a smooth function on $C=\mathbb{T}\times ]0,\infty [\times ]0,T[$ which decays rapidly when $y_{P}\rightarrow \infty$ . The indices $E$ , $P$ and $R$ denote, respectively, the Euler, Prandtl and remainder terms, and $y_{P}=\unicode[STIX]{x1D708}^{-1/2}y$ is the Prandtl variable. Note that, for simplicity, we have assumed that the flow is symmetric around the channel axis, so that the two $\unicode[STIX]{x1D714}_{P}$ terms correspond to two symmetric BLs of opposite sign at $y=0$ and $y=1$ .

By a classical multiple scales analysis, it can be formally shown that $\unicode[STIX]{x1D714}_{E}$ should satisfy the incompressible Euler equations in $\unicode[STIX]{x1D6FA}$ , and $\unicode[STIX]{x1D714}_{P}$ the Prandtl equations,

(2.16a ) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x2202}_{t}\unicode[STIX]{x1D714}_{P}+\unicode[STIX]{x1D735}(\boldsymbol{u}_{P}\unicode[STIX]{x1D714}_{P})=\unicode[STIX]{x2202}_{y_{P}}^{2}\unicode[STIX]{x1D714}_{P}, & \displaystyle\end{eqnarray}$$

(2.16b ) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D714}_{P}(x,y_{P},0)=0, & \displaystyle\end{eqnarray}$$

(2.16c ) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D713}_{P}(x,y_{P},t)=\int _{0}^{y_{P}}\text{d}y_{P}^{\prime }\int _{0}^{y_{P}^{\prime }}\text{d}y_{P}^{\prime \prime }\unicode[STIX]{x1D714}_{P}(x,y_{P}^{\prime \prime },t), & \displaystyle\end{eqnarray}$$

(2.16d ) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x2202}_{y_{P}}\unicode[STIX]{x1D714}_{P}(x,0,t)=-\unicode[STIX]{x2202}_{x}p_{E}(x,0,t), & \displaystyle\end{eqnarray}$$

$p_{E}$

$\unicode[STIX]{x1D714}_{E}$

$\unicode[STIX]{x1D714}$

(2.17) $$\begin{eqnarray}\displaystyle \int _{0}^{1}\widehat{\unicode[STIX]{x1D714}}_{Ek}(y)\exp (-|\unicode[STIX]{x1D705}|y)\,\text{d}y+\unicode[STIX]{x1D708}^{-1/2}\int _{0}^{1}\widehat{\unicode[STIX]{x1D714}}_{Pk}(y)\exp (-|\unicode[STIX]{x1D705}|y)\,\text{d}y=0, & & \displaystyle\end{eqnarray}$$

and by expressing the second integral with respect to $y_{P}$ and keeping the lowest-order term in $\unicode[STIX]{x1D708}$ , one has

(2.18) $$\begin{eqnarray}\displaystyle \int _{0}^{\infty }\widehat{\unicode[STIX]{x1D714}}_{Pk}(y_{P})\,\text{d}y_{P}=-\int _{0}^{1}\widehat{\unicode[STIX]{x1D714}}_{Ek}(y)\exp (-|\unicode[STIX]{x1D705}|y)\,\text{d}y. & & \displaystyle\end{eqnarray}$$

Then by integrating (2.16a ) over $[0,\infty ]$ , one finds that the contribution of the nonlinear term vanishes, and one is left with

(2.19) $$\begin{eqnarray}\displaystyle -\unicode[STIX]{x2202}_{t}\int _{0}^{1}\widehat{\unicode[STIX]{x1D714}}_{Ek}(y)\exp (-|\unicode[STIX]{x1D705}|y)\,\text{d}y=\int _{0}^{\infty }\unicode[STIX]{x2202}_{y_{P}}^{2}\widehat{\unicode[STIX]{x1D714}}_{Pk}(y)=-\unicode[STIX]{x2202}_{y_{P}}\widehat{\unicode[STIX]{x1D714}}_{Pk}(0), & & \displaystyle\end{eqnarray}$$

where, from the considerations in the preceding paragraph, it appears that the left-hand side can be identified with the pressure gradient $\unicode[STIX]{x2202}_{x}p_{E}(x,0,t)$ . Intuitively, the wall pressure gradient computed from the Euler solution creates vorticity at the boundary, which then diffuses inwards and evolves nonlinearly due to the flow it generates in the BL.

Since the Prandtl equations do not include diffusion parallel to the wall, in general nothing prevents the vorticity gradient in the $x$ direction from growing indefinitely – hence the possibility of finite-time singularity. More precisely, the mechanism proposed by Van Dommelen & Cowley (Reference Van Dommelen and Cowley1990) is that a fluid element is compressed to a point in the wall-parallel direction, and extends to infinity in the wall-normal direction. We shall denote by $t^{\ast }$ the time at which this first occurs, and by $x^{\ast }$ the corresponding $x$ location. From the scaling exponents computed by Van Dommelen & Cowley (Reference Van Dommelen and Cowley1990), we can deduce that, if the initial data are analytic, the spectrum of the solution will fill when approaching singularity with a characteristic cutoff parallel wavenumber scaling like

(2.20) $$\begin{eqnarray}\displaystyle k_{C}\propto (t^{\ast }-t)^{-3/2}. & & \displaystyle\end{eqnarray}$$

2.3 Interactive boundary layer model

As the singularity builds up in the Prandtl solution, the corresponding Navier–Stokes solution adopts a quite different behaviour. As first explained by Elliott et al. (Reference Elliott, Smith and Cowley1983), the first divergence between the two solutions occurs when the outer potential flow generated by the BL vorticity creates a pressure gradient perturbation of order 1 at the wall, which in turn impacts the inward diffusion of vorticity. This effect generically starts to take place when

(2.21) $$\begin{eqnarray}\displaystyle t^{\ast }-t=O(Re^{-2/11}). & & \displaystyle\end{eqnarray}$$

A rigorous asymptotic description of this new effect would require the modification of the vorticity ansatz (2.15) with new BLs, both in $x$ and in $y$ , coming into play. To avoid such complications, we follow Peridier et al. (Reference Peridier, Smith and Walker1991b ) and consider the finite-Reynolds-number description called the interactive boundary layer (IBL) model, which simply consists in modifying the Prandtl equations to include the new large-scale interaction, but without trying to rescale the solution a priori. Ansatz (2.15) therefore remains valid, except that $\unicode[STIX]{x1D714}_{P}$ is replaced by $\unicode[STIX]{x1D714}_{I}$ , the solution of the interactive equations, which we shall now derive.

Since we are working with the vorticity formulation, we are blind to potential flow perturbations, but their effect manifests itself through the integral boundary condition (2.12) on $\unicode[STIX]{x1D714}$ . Starting again from (2.17), but expanding the exponential up to order $Re^{-1/2}$ , yields

(2.22) $$\begin{eqnarray}\displaystyle \int _{0}^{1}\widehat{\unicode[STIX]{x1D714}}_{Ik}(y_{P})\,\text{d}y_{P}-\unicode[STIX]{x1D708}^{1/2}|\unicode[STIX]{x1D705}|\int _{0}^{\infty }y_{P}\widehat{\unicode[STIX]{x1D714}}_{Ik}(y_{P})\,\text{d}y_{P}=-\int _{0}^{1}\widehat{\unicode[STIX]{x1D714}}_{Ek}(y)\exp (-|\unicode[STIX]{x1D705}|y)\,\text{d}y, & & \displaystyle\end{eqnarray}$$

and, following the same procedure as above, leads to a perturbed boundary condition for $\unicode[STIX]{x1D714}_{I}$ ,

(2.23) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x2202}_{y_{P}}\widehat{\unicode[STIX]{x1D714}}_{Ik}(0)=-\text{i}\unicode[STIX]{x1D705}\widehat{p}_{Ek}(0)-\unicode[STIX]{x1D708}^{1/2}|\unicode[STIX]{x1D705}|\unicode[STIX]{x2202}_{t}\widehat{\unicode[STIX]{x1D6FD}}_{Ik}, & & \displaystyle\end{eqnarray}$$

where

(2.24) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6FD}_{Ik}(x)=\int _{0}^{\infty }y_{P}\unicode[STIX]{x1D714}_{I}(x,y_{P})\,\text{d}y_{P}. & & \displaystyle\end{eqnarray}$$

On the other hand, by multiplying (2.16a ) by $y_{P}$ and integrating over $y_{P}$ , we obtain an expression for $\unicode[STIX]{x2202}_{t}\unicode[STIX]{x1D6FD}_{I}$ , which closes the problem.

As a side remark, let us note that the classical name ‘interactive boundary layer’ for this model is misleading, since in fact no retroaction of the Prandtl layer onto the bulk Euler flow is taken into account. An alternative name could be ‘wet boundary layer’, which better encompasses the notion that the potential far flow affects the boundary layer equations only through a passive effect.

2.4 Orr–Sommerfeld model

Several numerical studies suggest that a linear instability mechanism could play a role in the detachment process. Since we are concerned with an unsteady problem, the notion of linear instability should be understood here in an asymptotic sense, in terms of a rescaled time variable in which the evolution of the base flow can be neglected. Moreover, since we are looking for an instability happening at high wavenumbers in the parallel direction, we also neglect, for the time being, the parallel variation of the base flow, or in other words we study the possible occurrence of perturbations which have a large parallel wavenumber $k$ compared to the characteristic parallel wavenumber $L^{-1}$ of the flow prior to detachment. The combination of both hypotheses constitutes the frozen flow approximation. Its domain of validity could be properly evaluated only by resorting to a multiple-time-scale asymptotic analysis, which we have not yet achieved in this setting.

2.4.1 Formulation

Under these two simplifying hypotheses, we are brought back to Rayleigh’s classical shear flow stability problem, later generalized to viscous fluids by Orr and Sommerfeld. In the case of a boundary layer, several simplifications are possible which allowed Tollmien (Reference Tollmien1929) to obtain an elegant asymptotic description of the modes now known as Tollmien–Schlichting waves, and of the corresponding stability region in the $(Re,k)$ plane, which was later confirmed experimentally by Schubauer & Skramstad (Reference Schubauer and Skramstad1947). For a more recent review on the subject, see Reed, Saric & Arnal (Reference Reed, Saric and Arnal1996). Although most of the material presented here is classical (see Lin Reference Lin1967), previous studies have mostly emphasized the computation of the critical Reynolds number, so that it is instructive to rederive the main results directly in the $Re\rightarrow \infty$ limit, which concerns us here.

For small perturbations $\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D713}(x,y_{P},t_{2})=\unicode[STIX]{x1D719}(y_{P})\text{e}^{\text{i}(\unicode[STIX]{x1D705}x-\unicode[STIX]{x1D6FC}t_{2})}$ to the streamfunction, the profile function $\unicode[STIX]{x1D719}$ satisfies the Orr–Sommerfeld equation

(2.25) $$\begin{eqnarray}\displaystyle (u_{P}-c)(\unicode[STIX]{x1D719}^{\prime \prime }-\unicode[STIX]{x1D708}\unicode[STIX]{x1D705}^{2}\unicode[STIX]{x1D719})-u_{P}^{\prime \prime }\unicode[STIX]{x1D719}=\frac{1}{\text{i}\unicode[STIX]{x1D705}}(\unicode[STIX]{x1D719}^{\prime \prime \prime \prime }-\unicode[STIX]{x1D708}\unicode[STIX]{x1D705}^{2}\unicode[STIX]{x1D719}^{\prime \prime }+\unicode[STIX]{x1D708}^{2}\unicode[STIX]{x1D705}^{4}\unicode[STIX]{x1D719}), & & \displaystyle\end{eqnarray}$$

where $c=\unicode[STIX]{x1D6FC}/\unicode[STIX]{x1D705}$ is the phase velocity, and primes denote derivatives with respect to $y_{P}$ . Note that $c$ is in general a complex number, and unstable perturbations are those for which $c$ has a strictly positive imaginary part. Now assuming that

(2.26) $$\begin{eqnarray}\displaystyle L^{-1}\ll k\ll L^{-1}Re^{1/2}, & & \displaystyle\end{eqnarray}$$

(2.25) simplifies to

(2.27) $$\begin{eqnarray}\displaystyle (u_{P}-c)\unicode[STIX]{x1D719}^{\prime \prime }-u_{P}^{\prime \prime }\unicode[STIX]{x1D719}=\frac{1}{\text{i}\unicode[STIX]{x1D705}}\unicode[STIX]{x1D719}^{\prime \prime \prime \prime }. & & \displaystyle\end{eqnarray}$$

The no-slip boundary condition translates to $\unicode[STIX]{x1D719}(0)=\unicode[STIX]{x1D719}^{\prime }(0)=0$ . Assuming from now on that $k>0$ without loss of generality, the boundary condition for $y_{P}\rightarrow +\infty$ can be obtained by matching $\unicode[STIX]{x1D719}$ with a harmonic outer solution of the form $\exp (-\unicode[STIX]{x1D705}y)$ using the hypothesis that vorticity vanishes outside the BL, which means that

(2.28) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D719}(y_{P})\underset{y_{P}\rightarrow \infty }{=}A(1-\unicode[STIX]{x1D705}\unicode[STIX]{x1D708}^{1/2}y_{P})+o(\unicode[STIX]{x1D705}\unicode[STIX]{x1D708}^{1/2}y_{P}). & & \displaystyle\end{eqnarray}$$

Note that it is essential to keep the first-order term in this expression in order to find unstable modes. Following Tollmien (Reference Tollmien1929), we now deal with the singular perturbation problem (2.27) by first considering inviscid solutions, and then adding a boundary layer.

2.4.2 Inviscid mode

Neglecting the viscous contribution in (2.27), we obtain

(2.29) $$\begin{eqnarray}\displaystyle (u_{P}-c)\unicode[STIX]{x1D719}^{\prime \prime }-u_{P}^{\prime \prime }\unicode[STIX]{x1D719}=0, & & \displaystyle\end{eqnarray}$$

which admits the obvious regular solution

(2.30) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D719}_{1,c}=u_{P}-c, & & \displaystyle\end{eqnarray}$$

but is singular at any point where $u_{P}=c$ . To construct another independent solution $\unicode[STIX]{x1D719}_{2,c}$ , we now assume that $c$ does not lie directly on the real axis, and we make the change of unknown $\unicode[STIX]{x1D719}=(u_{P}-c)f,$ leading to

(2.31) $$\begin{eqnarray}\displaystyle (u_{P}-c)f^{\prime \prime }+2u_{P}^{\prime }f^{\prime }=0 & & \displaystyle\end{eqnarray}$$

and thus

(2.32) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D719}_{2,c}(y_{P})=(u_{P}-c)\left(\int _{\infty }^{y_{P}}\left(\left(\frac{u_{\infty }-c}{u_{P}-c}\right)^{2}-1\right)+y_{P}\right), & & \displaystyle\end{eqnarray}$$

where $u_{\infty }$ is the velocity outside the boundary layer. By combining $\unicode[STIX]{x1D719}_{1,c}$ and $\unicode[STIX]{x1D719}_{2,c}$ , a solution $\unicode[STIX]{x1D719}_{out}$ satisfying the condition (2.28) at $+\infty$ is readily obtained:

(2.33) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D719}_{out}=\unicode[STIX]{x1D719}_{1,c}-\unicode[STIX]{x1D705}\unicode[STIX]{x1D708}^{1/2}\unicode[STIX]{x1D719}_{2,c}. & & \displaystyle\end{eqnarray}$$

2.4.3 Viscous correction

We now look for a viscous sublayer correction $\unicode[STIX]{x1D719}_{in}$ which is necessary since $\unicode[STIX]{x1D719}_{out}$ does not in general satisfy the no-slip boundary condition at $y_{P}=0$ . For small $y_{P}$ , (2.27) reduces to

(2.34) $$\begin{eqnarray}\displaystyle u_{P}^{\prime }(z_{0}(c))(y_{P}-z_{0}(c))\unicode[STIX]{x1D719}^{\prime \prime }-u_{P}^{\prime \prime }(0)\unicode[STIX]{x1D719}(0)=\frac{1}{\text{i}\unicode[STIX]{x1D705}}\unicode[STIX]{x1D719}^{\prime \prime \prime \prime }, & & \displaystyle\end{eqnarray}$$

where we have defined $z_{0}(c)$ to be the solution with the smallest real part to the equation $u(z)=c$ (see appendix A). An inner variable can then be defined in the viscous sublayer by $\unicode[STIX]{x1D702}=\unicode[STIX]{x1D700}(y_{P}-z(c))|\unicode[STIX]{x1D705}u_{P}^{\prime }(z(c))|5^{1/3}$ , where $\unicode[STIX]{x1D700}=\text{sign}(u_{P}^{\prime }(z(c)))$ , leading to, with $\unicode[STIX]{x1D705}\gg 1$ ,

(2.35) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D702}\widetilde{\unicode[STIX]{x1D719}}^{\prime \prime }=-\text{i}\widetilde{\unicode[STIX]{x1D719}}^{\prime \prime \prime \prime }. & & \displaystyle\end{eqnarray}$$

We are interested in a solution of this equation which remains bounded and whose derivative tends to zero when $y_{P}\rightarrow \infty$ . When $\unicode[STIX]{x1D700}>0$ , this limit is equivalent to $\unicode[STIX]{x1D702}\rightarrow \infty$ , and a solution to the problem was given by Tollmien (Reference Tollmien1929) in terms of Hankel functions:

(2.36) $$\begin{eqnarray}\displaystyle \widetilde{\unicode[STIX]{x1D719}}_{in}(\unicode[STIX]{x1D702})=\int _{\unicode[STIX]{x1D702}}^{\infty }\text{d}\unicode[STIX]{x1D702}^{\prime }\int _{\unicode[STIX]{x1D702}^{\prime }}^{\infty }\text{d}\unicode[STIX]{x1D702}^{\prime \prime }\,{\unicode[STIX]{x1D702}^{\prime \prime }}^{1/2}\text{H}_{1/3}^{(1)}({\textstyle \frac{2}{3}}(\text{i}\unicode[STIX]{x1D702}^{\prime \prime })^{3/2}). & & \displaystyle\end{eqnarray}$$

Note that, as long as it is expressed in terms of the $\unicode[STIX]{x1D702}$ variable, this solution is universal. Expressed as a function of $y_{P}$ , it reads

(2.37) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D719}_{in}(y_{P})=\widetilde{\unicode[STIX]{x1D719}}_{in}\left(\left(y_{P}-\frac{c}{u_{P}^{\prime }(0)}\right)|\unicode[STIX]{x1D705}u_{P}^{\prime }(0)|^{1/3}\right). & & \displaystyle\end{eqnarray}$$

In the case $\unicode[STIX]{x1D700}<0$ , $y_{P}\rightarrow \infty$ corresponds to $\unicode[STIX]{x1D702}\rightarrow -\infty$ , and the solution $\widetilde{\unicode[STIX]{x1D719}}_{in}(\unicode[STIX]{x1D702})$ should be adapted accordingly.

2.4.4 Construction of unstable modes

Now, in the asymptotic regime, any admissible solution $\unicode[STIX]{x1D719}$ of (2.27) can approximately be expressed as

(2.38) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D719}=A\unicode[STIX]{x1D719}_{out}+B\unicode[STIX]{x1D719}_{in}, & & \displaystyle\end{eqnarray}$$

so that the boundary conditions $\unicode[STIX]{x1D719}(0)=\unicode[STIX]{x1D719}^{\prime }(0)=0$ translate to the linear system

(2.39) $$\begin{eqnarray}\displaystyle \left\{\begin{array}{@{}l@{}}A\unicode[STIX]{x1D719}_{out}(0)+B\unicode[STIX]{x1D719}_{in}(0)=0,\quad \\ A\unicode[STIX]{x1D719}_{out}^{\prime }(0)+B\unicode[STIX]{x1D719}_{in}^{\prime }(0)=0.\quad \end{array}\right. & & \displaystyle\end{eqnarray}$$

In order for a non-trivial solution to exist, this system should be degenerate, i.e.

(2.40) $$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x1D719}_{in}(0)}{\unicode[STIX]{x1D719}_{in}^{\prime }(0)}=\frac{\unicode[STIX]{x1D719}_{out}(0)}{\unicode[STIX]{x1D719}_{out}^{\prime }(0)}. & & \displaystyle\end{eqnarray}$$

On the one hand, denoting

(2.41) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D702}_{w}=-c\unicode[STIX]{x1D705}^{1/3}|u_{P}^{\prime }(0)|^{-2/3}, & & \displaystyle\end{eqnarray}$$

the position of the wall in terms of the $\unicode[STIX]{x1D702}$ variable, it is shown following Tollmien (Reference Tollmien1929) that

(2.42) $$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x1D719}_{in}(0)}{\unicode[STIX]{x1D719}_{in}^{\prime }(0)}=-\frac{c}{u_{P}^{\prime }(0)}F(-\unicode[STIX]{x1D702}_{w})\quad \text{if }u_{P}^{\prime }(0)>0, & & \displaystyle\end{eqnarray}$$

where $F$ is a one-parameter complex function known as the Tietjens function $F$ . Although there is no closed analytical formula for $F$ , it is easily approximated from (2.37) using quadrature formulas. A graphical representation is given in figure 1.

Figure 1.

Imaginary versus real part of the Tietjens function $F(\unicode[STIX]{x1D702})$ for $\unicode[STIX]{x1D702}\in [2,20]$ . The sampling points (every $0.1$ ) are indicated by crosses.

On the other hand, denoting for convenience of notation

(2.43) $$\begin{eqnarray}\displaystyle I_{c}=-c\int _{0}^{\infty }\left(\frac{1}{(u_{P}(y)-c)^{2}}-\frac{1}{(u_{\infty }-c)^{2}}\right)\,\text{d}y & & \displaystyle\end{eqnarray}$$

and using (2.33), it is shown that

(2.44) $$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x1D719}_{out}(0)}{\unicode[STIX]{x1D719}_{out}^{\prime }(0)}=-\frac{c}{u_{P}^{\prime }(0)}E(\unicode[STIX]{x1D705}\unicode[STIX]{x1D708}^{1/2},c), & & \displaystyle\end{eqnarray}$$

where

(2.45) $$\begin{eqnarray}\displaystyle E(\unicode[STIX]{x1D705}\unicode[STIX]{x1D708}^{1/2},c)=1-\frac{\unicode[STIX]{x1D705}\unicode[STIX]{x1D708}^{1/2}(u_{\infty }-c)^{2}}{cu_{P}^{\prime }(0)\left(1-\displaystyle \frac{\unicode[STIX]{x1D705}\unicode[STIX]{x1D708}^{1/2}(u_{\infty }-c)^{2}}{c}\left(I_{c}-\displaystyle \frac{1}{u^{\prime }(0)}\right)\right)}. & & \displaystyle\end{eqnarray}$$

Injecting (2.42) and (2.44) back into the degeneracy condition (2.40), we obtain the dispersion relation

(2.46) $$\begin{eqnarray}\displaystyle F(-\unicode[STIX]{x1D702}_{w})=E(\unicode[STIX]{x1D705}\unicode[STIX]{x1D708}^{1/2},c). & & \displaystyle\end{eqnarray}$$

The profile is linearly unstable if and only if this equation has solutions such that $\text{Im}(c)>0$ . In the range of $\unicode[STIX]{x1D705}$ we are considering, there can be no solutions if $c$ is of order $u_{\infty }$ or larger, because then $\unicode[STIX]{x1D702}_{w}\rightarrow -\infty$ , and therefore $F(-\unicode[STIX]{x1D702}_{w})\rightarrow 0$ , whereas $E$ remains of order $1$ . In the following, we therefore look for solutions under the restriction that $c\ll u_{\infty }$ .

The asymptotic behaviour of $I_{c}$ for small $c$ is dominated by what happens near solutions of $u(y_{P})=0$ . As established in appendix A, if $u_{P}$ and $u_{P}^{\prime }$ do not have any zeros in common, then

(2.47) $$\begin{eqnarray}\displaystyle I_{c}\underset{c\rightarrow 0}{=}\frac{1}{u_{P}^{\prime }(0)}-c\frac{u_{P}^{\prime \prime }(0)}{u_{P}^{\prime }(0)^{3}}\ln \left(-\frac{c}{u_{P}^{\prime }(0)}\right)-c\text{i}\unicode[STIX]{x03C0}\unicode[STIX]{x1D6E5}_{u}+cK_{R}+o(c), & & \displaystyle\end{eqnarray}$$

where $K_{R}$ is a real constant, $\ln$ is the principal branch of the complex logarithm, and

(2.48) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6E5}_{u}=\mathop{\sum }_{\{y_{P}>0\mid u_{P}(y_{P})=0\}}\frac{u_{P}^{\prime \prime }(y_{P})}{u_{P}^{\prime }(y_{P})^{3}}. & & \displaystyle\end{eqnarray}$$

The behaviour of this asymptotic expression when $c$ approaches the real axis is tricky and should be considered with care. If $\text{Re}(c/u_{P}^{\prime }(0))<0$ , the argument of the complex logarithm lies on the right-hand side of the complex plane, and the limit $\text{Im}(c)\rightarrow 0$ is well behaved. But if $\text{Re}(c/u_{P}^{\prime }(0))>0$ , the limit does not exist strictly speaking, and $I_{c}$ should be considered as multi-valued, which is well captured by replacing $\unicode[STIX]{x1D6E5}_{u}$ by

(2.49) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6E5}_{u}^{\pm }=\left\{\begin{array}{@{}ll@{}}\unicode[STIX]{x1D6E5}_{u}\quad & \text{ if }\text{Re}(c/u_{P}^{\prime }(0))<0,\\ \displaystyle \pm \frac{u_{P}^{\prime \prime }(0)}{u_{P}^{\prime }(0)^{3}}+\unicode[STIX]{x1D6E5}_{u}\quad & \text{ if }\text{Re}(c/u_{P}^{\prime }(0))>0.\end{array}\right. & & \displaystyle\end{eqnarray}$$

By inserting our estimate for $I_{c}$ into (2.45), and assuming for simplicity that $\unicode[STIX]{x1D705}\unicode[STIX]{x1D708}^{1/2}\log |c|\ll 1$ , we obtain the following estimates for the real and imaginary parts of $E$ (respectively) when $c$ is close to the real axis:

(2.50a ) $$\begin{eqnarray}\displaystyle & \displaystyle \text{Re}(E)\sim 1-\frac{\unicode[STIX]{x1D705}\unicode[STIX]{x1D708}^{1/2}u_{\infty }^{2}}{cu_{P}^{\prime }(0)}, & \displaystyle\end{eqnarray}$$

(2.50b ) $$\begin{eqnarray}\displaystyle & \displaystyle \text{Im}(E)\sim \frac{\unicode[STIX]{x1D705}^{2}\unicode[STIX]{x1D708}u_{\infty }^{4}\unicode[STIX]{x1D6E5}_{u}^{\pm }\unicode[STIX]{x03C0}}{cu_{P}^{\prime }(0)}, & \displaystyle\end{eqnarray}$$

where we have assumed that $\unicode[STIX]{x1D6E5}_{u}^{\pm }\neq 0$ (this analysis therefore does not apply directly to the Blasius boundary layer). Since the imaginary part of $E$ is negligible compared to its real part, (2.40) can be satisfied only in the neighbourhood of points where $F$ is purely real. This happens for $\unicode[STIX]{x1D702}_{w}\rightarrow -\infty$ , as well as for a certain finite value $\unicode[STIX]{x1D702}_{0}\simeq -2.3$ where $F(-\unicode[STIX]{x1D702}_{0})\simeq 0.56$ (see figure 1).

In the case $u_{P}^{\prime }(0)<0$ , the latter leads, with (2.50a ), to a single critical wavenumber

(2.51) $$\begin{eqnarray}\displaystyle k_{1}\simeq 0.161|u_{P}^{\prime }(0)|^{5/4}u_{\infty }^{-3/2}\unicode[STIX]{x1D708}^{-3/8}, & & \displaystyle\end{eqnarray}$$

beyond which the modes are unstable. In the case $u_{P}^{\prime }(0)>0$ , the multi-valuedness of $I_{c}$ implies that $k_{1}$ splits into two critical wavenumbers with very close values, and, therefore, a negligible instability region.

The critical point near $\unicode[STIX]{x1D702}_{w}\rightarrow -\infty$ corresponds to $F(-\unicode[STIX]{x1D702}_{w})\rightarrow 0$ , so from (2.50a )

(2.52) $$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x1D705}\unicode[STIX]{x1D708}^{1/2}u_{\infty }^{2}}{cu_{P}^{\prime }(0)}\simeq 1. & & \displaystyle\end{eqnarray}$$

To obtain another equation relating $c$ and $\unicode[STIX]{x1D705}$ , we use the estimate $F(-\unicode[STIX]{x1D702}_{w}){\sim}_{\unicode[STIX]{x1D702}\rightarrow -\infty }-\text{e}^{\text{i}\unicode[STIX]{x03C0}/4}|\unicode[STIX]{x1D702}|^{-3/2},$ which, combined with (2.50b ), leads to

(2.53) $$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x1D705}^{2}\unicode[STIX]{x1D708}u_{\infty }^{4}\unicode[STIX]{x1D6E5}_{u}^{\pm }\unicode[STIX]{x03C0}}{cu_{P}^{\prime }(0)}\simeq \frac{\sqrt{2}}{2}|\unicode[STIX]{x1D702}_{w}|^{-3/2}. & & \displaystyle\end{eqnarray}$$

This equation has a simple root if and only if $\unicode[STIX]{x1D6E5}_{u}>0$ and $u_{P}^{\prime }(0)<0$ , in which case we obtain a second critical wavenumber

(2.54) $$\begin{eqnarray}\displaystyle k_{2}=0.0968|u_{P}^{\prime }(0)|^{1/2}|\unicode[STIX]{x1D6E5}_{u}|^{-1/3}|u_{\infty }|^{-5/3}\unicode[STIX]{x1D708}^{-5/12}. & & \displaystyle\end{eqnarray}$$

In other cases, if there exists an unstable region for $k>k_{1}$ , its upper bound cannot be found under our current restriction $k\ll \unicode[STIX]{x1D708}^{-1/2}$ , which implies that it extends at least up to wavenumbers scaling like $\unicode[STIX]{x1D708}^{-1/2}$ , corresponding to what is usually called the Rayleigh instability. This observation should be kept in mind as it is one of the key elements of the detachment scenario we will propose below.

Another important quantity we need to estimate is the growth rate of the unstable modes. From (2.45), we see that, since $c$ remains small in absolute value, the growth rate $\unicode[STIX]{x1D6FC}$ of the instability satisfies

(2.55) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6FC}\sim \frac{u_{\infty }^{2}\text{Im}(F(-\unicode[STIX]{x1D702}_{w}))}{|u_{P}^{\prime }(0)|}\unicode[STIX]{x1D708}^{1/2}\unicode[STIX]{x1D705}^{2}. & & \displaystyle\end{eqnarray}$$

To sum up, the generic instability expected to play a role in such boundary layer flows in the inviscid limit manifests itself by the growth of wavepackets in the vicinity of the boundary confined in physical space to regions where $u_{P}^{\prime }(0)<0$ (recirculation bubbles), and whose parallel wavenumber support extends from $O(Re^{3/8})$ to at least $O(Re^{1/2})$ .

2.4.5 The case $u_{P}^{\prime }(0)=0$

To be complete, our analysis should also take into account the case $u_{P}^{\prime }(0)=0$ , investigated in detail by Hughes & Reid (Reference Hughes and Reid1965). Going back to the general expression (2.33) for the outer solution, we obtain in the case $u_{P}^{\prime }(0)=0$ that

(2.56) $$\begin{eqnarray}\displaystyle E(\unicode[STIX]{x1D705}\unicode[STIX]{x1D708}^{1/2},c)=\frac{\unicode[STIX]{x1D719}_{out}(0)}{\unicode[STIX]{x1D719}_{out}^{\prime }(0)}=-\frac{c^{2}}{\unicode[STIX]{x1D705}\unicode[STIX]{x1D708}^{1/2}(u_{\infty }-c)^{2}}-cI_{c} & & \displaystyle\end{eqnarray}$$

or, with (A 18), and when $c$ is close to the real axis,

(2.57a ) $$\begin{eqnarray}\displaystyle & \displaystyle \text{Re}(E)\simeq -\frac{c^{2}}{\unicode[STIX]{x1D705}\unicode[STIX]{x1D708}^{1/2}u_{\infty }^{2}}, & \displaystyle\end{eqnarray}$$

(2.57b ) $$\begin{eqnarray}\displaystyle & \displaystyle \text{Im}(E)\simeq -\frac{\unicode[STIX]{x03C0}}{4}\left(\frac{2c}{u_{P}^{\prime \prime }(0)}\right)^{1/2}. & \displaystyle\end{eqnarray}$$

Concerning the inner solution, (2.42) is replaced by

(2.58) $$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x1D719}_{in}(0)}{\unicode[STIX]{x1D719}_{in}^{\prime }(0)}=-z(c)F(-\unicode[STIX]{x1D702}_{w}), & & \displaystyle\end{eqnarray}$$

which, combined with (2.57b ), yields

(2.59) $$\begin{eqnarray}\displaystyle \text{Im}(F(-\unicode[STIX]{x1D702}_{w}))=-\frac{\unicode[STIX]{x03C0}}{4}, & & \displaystyle\end{eqnarray}$$

or equivalently, according to Hughes & Reid (Reference Hughes and Reid1965),

(2.60a,b ) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D702}_{w}\simeq -0.488,\quad \text{Re}(F(-\unicode[STIX]{x1D702}_{w}))\simeq 1.580, & & \displaystyle\end{eqnarray}$$

and therefore

(2.61) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{c^{3/2}}{\unicode[STIX]{x1D705}\unicode[STIX]{x1D708}^{1/2}u_{\infty }^{2}}=\sqrt{\frac{2}{u_{P}^{\prime \prime }(0)}}1.580, & \displaystyle\end{eqnarray}$$

(2.62) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{2c}{u_{P}^{\prime \prime }(0)}\unicode[STIX]{x1D705}^{1/3}(2cu^{\prime \prime }(0))^{1/3}=0.488, & \displaystyle\end{eqnarray}$$

which finally gives us the critical wavenumber

(2.63) $$\begin{eqnarray}\displaystyle k_{1}=0.0279u_{P}^{\prime \prime }(0)^{10/11}u_{\infty }^{16/11}\unicode[STIX]{x1D708}^{-4/11}. & & \displaystyle\end{eqnarray}$$

The growth rate, obtained following the same reasoning as above with (2.45) replaced by (2.56), reads

(2.64) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6FC}\sim u_{\infty }\text{Im}(F(-\unicode[STIX]{x1D702}_{w}))\unicode[STIX]{x1D708}^{1/4}\unicode[STIX]{x1D705}^{3/2}. & & \displaystyle\end{eqnarray}$$

2.4.6 Physical interpretation

In this section we formulate some conjectures regarding the physical interpretation of the above model. We have shown that, subject to the validity of the frozen flow approximation, all BL flows containing recirculation bubbles are subject to Tollmien–Schlichting–Rayleigh instabilities for wavenumbers $k\in [k_{1},k_{2}]$ , where $k_{1}$ and $k_{2}$ both diverge to $\infty$ when $Re\rightarrow \infty$ .

Therefore a plausible scenario for detachment may begin as follows. Suppose that initially the flow is very smooth, for example, that it has analytic regularity, i.e. its Fourier coefficients decay exponentially with $k$ , and that a recirculation bubble appears due to the Prandtl BL nonlinear dynamics. Our analysis then suggests that within the recirculation bubble the range $[k_{1},k_{2}]$ is subject to a fast linear instability. Note that, since we have found that $k_{1}$ , $k_{2}$ and the growth rate $\unicode[STIX]{x1D6FC}$ (see (2.55)) diverge with $Re$ , there are reasonable chances that the frozen flow approximation becomes more and more accurate for higher $Re$ . However, for sufficiently large Reynolds number, the initial excitation of such high-wavenumber modes is so small that they do not have time to grow and the Prandtl solution remains a good approximation.

But if and when a Prandtl singularity builds up, it starts feeding non-negligible excitations into the interval $[k_{1},k_{2}]$ . In the competition between the oncoming singularity and the growth of unstable modes, it is interesting to determine which modes first reach a finite amplitude, and when this occurs.

Now if we replace $\unicode[STIX]{x1D705}$ by the characteristic excitation (2.20) generated by the Prandtl dynamics some time $t^{\ast }-t$ before the singularity, we obtain

(2.65) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6FC}\propto \unicode[STIX]{x1D708}^{1/2}(t^{\ast }-t)^{-3}. & & \displaystyle\end{eqnarray}$$

With this growth rate, the first perturbations to reach order 1 occur at a time

(2.66) $$\begin{eqnarray}\displaystyle t^{\ast }-t=O(Re^{-1/4}). & & \displaystyle\end{eqnarray}$$

By comparing this result with (2.21), we note that this occurs later than the perturbations due to large-scale interactions, as described by the IBL model. Therefore, the BL profile resulting from an IBL computation, not the Prandtl profile, should be used as base profile when performing the stability analysis. This confirms the analysis of Gargano et al. (Reference Gargano, Sammartino, Sciacca and Cassel2014), who pointed out that what they call a large-scale interaction always precedes the approach to detachment.

In the region with reversed flow near the wall, the width of the unstable wavenumber range scales likes $O(Re^{1/2})$ . Assuming that all the modes grow simultaneously and reach order 1, this means that the support of $\widehat{\unicode[STIX]{x1D714}}$ extends to $k\propto O(Re^{1/2})$ , while the amplitude of the modes continues to scale like $O(Re^{1/2})$ . Owing to the properties of the inverse Fourier transform, these scalings immediately imply that the profile of $\unicode[STIX]{x1D714}$ very near the wall has a kind of wavepacket shape with amplitude scaling like $O(Re^{1})$ .

During the linear phase, the characteristic wall-normal extent of such modes is controlled by the considerations of § 2.4.3, i.e. $\unicode[STIX]{x1D705}^{-1/3}Re^{-1/2}\sim Re^{-2/3}$ . But once the unstable modes have reached order 1 and the amplitude of $\unicode[STIX]{x1D714}$ scales like $O(Re^{1})$ (due to the superposition of all modes as noted above), nonlinear vorticity advection effects imply that the characteristic scale becomes $O(Re^{-1})$ , which gives us a possible physical explanation for the Kato layer.

3.1 Set-up

To trigger an unsteady separation process, we have chosen an initial configuration inspired by the dipole of Orlandi (Reference Orlandi1990), later modified by Clercx & van Heijst (Reference Clercx and van Heijst2002). However, this dipole has the drawback of generating a secondary, weaker dipole propagating in the opposite direction which is computed at a waste. For efficiency reasons, we have thus preferred a quadrupole configuration, which is symmetric both around the channel axis and around the midplane, thus sparing three-quarters of the domain size for a given $Re$ . It is defined in terms of its streamfunction as follows:

(3.1) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D713}_{i}(x,y)=Axy\exp \left(-\frac{(x-x_{0})^{2}+(y-y_{0})^{2}}{2s^{2}}\right), & & \displaystyle\end{eqnarray}$$

where $A=0.625847306637464$ and $s=6.3661977236758$ determine the amplitude of the vortices and their size, and $(x_{0},y_{0})=(0.5,0.5)$ their initial location. Note that the boundary conditions are satisfied only approximately by this velocity field, but in fact

(3.2) $$\begin{eqnarray}\displaystyle v_{i}(x,y=0)\approx 10^{-15}, & & \displaystyle\end{eqnarray}$$

which is in any case of the same order as the round-off error in double-precision arithmetic.

Owing to the symmetry of this initial condition, the analysis can be restricted without loss of information to the subdomain $K=[0,1/2]\times [0,1/2]$ . The streamlines of $\boldsymbol{u}_{i}$ in $K$ are shown in figure 2. The definitions and initial values of several integral quantities which will be important in our study are given in table 1.

Figure 2.

Streamlines of initial velocity field in the subdomain $K=[1/2,1]\times [0,1/2]$ .

Table 1.

Definitions and initial values of several integral quantities of interest.

In this work we shall analyse the flows obtained by solving the Navier–Stokes equations numerically up to $T=57.05$ for $\unicode[STIX]{x1D708}$ decreasing from $4\times 10^{-7}$ to $5\times 10^{-8}$ by factors of $\sqrt{2}$ (i.e. seven different values in total). Reynolds numbers corresponding to these values of $\unicode[STIX]{x1D708}$ are defined according to (2.2), where $U\simeq 2.42\times 10^{-2}$ is the initial maximum velocity, and $L=2s\simeq 1.27\times 10^{-1}$ is a measure of the size of the quadrupole. Both $\unicode[STIX]{x1D708}$ and $Re$ are provided up to three significant digits in table 2. To facilitate comparison with previous results concerning dipole–wall collisions (see e.g. Clercx & van Heijst Reference Clercx and van Heijst2017), we have also included the Reynolds number $Re_{rms}$ computed from the root-mean-square (r.m.s.) velocity $U_{rms}=4.25\times 10^{-3}$ and channel width instead, which is of the same order of magnitude and a factor 1.4 larger.

Table 2.

Parameters of numerical experiments. All figures in this table are given up to three significant digits.

3.2 Navier–Stokes solver

To solve the initial value problem for the Navier–Stokes equations, derivatives in the periodic $x$ direction are computed with spectral resolution from their sine and cosine series expansions. As we shall see below, we will need two different grid refinements in the $x$ direction, with $N_{x,1}$ (respectively $N_{x,2}$ ) collocation points and a corresponding resolution of $\unicode[STIX]{x1D6E5}_{x,1}$ (respectively $\unicode[STIX]{x1D6E5}_{x,2}$ ). Since the $y_{P}$ direction is not periodic, derivatives in the $y_{P}$ direction have to be treated differently. The Chebyshev scheme in the wall-normal direction is accurate but very costly and requires Gauss collocation points, which are optimal from an approximation theory point of view. However, according to our linear analysis, Gauss collocation points are not optimal. Indeed, to have sufficient resolution in the bulk flow, the number of modes should scale at least like $Re^{1/2}$ , which in turn imposes a wall-normal resolution $O(Re^{-1})$ in the whole Prandtl boundary layer. But this would be a waste of resolution since the initial growth of a wavepacket in the event of a Tollmien–Schlichting instability would occur only in a sublayer of thickness $Re^{-2/3}$ . We have therefore preferred to turn to fifth-order compact finite differences (Lele Reference Lele1992; Gamet et al. Reference Gamet, Ducros, Nicoud and Poinsot1999). Denoting by $f_{i}$ approximate values of a function $f$ on the uniform grid defined by $y_{1,i}=i/(N_{y_{P}}-1)~(0\leqslant i<N_{y_{P}})$ , and by $f_{i}^{\prime }$ and $f_{i}^{\prime \prime }$ approximations of its first and second derivatives at the same locations, we impose fifth-order accuracy by requiring that

(3.3a ) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D6FC}_{i}f_{i-1}^{\prime }+f_{i}^{\prime }+\unicode[STIX]{x1D6FD}_{i}f_{i+1}^{\prime }=A_{i}f_{i-1}+B_{i}f_{i}+C_{i}f_{i+1}, & \displaystyle\end{eqnarray}$$

(3.3b ) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D6FE}_{i}f_{i-1}^{\prime \prime }+f_{i}^{\prime \prime }+\unicode[STIX]{x1D6FF}_{i}f_{i+1}^{\prime \prime }=D_{i}f_{i-1}+E_{i}f_{i}+F_{i}f_{i+1}, & \displaystyle\end{eqnarray}$$

Note, however, that these expressions are only valid for $1\leqslant i<N_{y}-1$ , so that two additional equations are needed to determine $f_{i}^{\prime }$ and $f_{i}^{\prime \prime }$ uniquely. For the computation of $\unicode[STIX]{x2202}_{y}(v\unicode[STIX]{x1D714})$ , they are obtained by noting that the derivative vanishes at $y=0$ and $y=1$ , which is a direct consequence of the boundary conditions on $u$ and $v$ and of incompressibility.

For the viscous term $\unicode[STIX]{x2202}_{y}^{2}\unicode[STIX]{x1D714}$ , they should follow from the boundary conditions (2.14). To derive them, the integrals are first discretized by a fifth-order local quadrature formula. To preserve accuracy, $\unicode[STIX]{x1D714}$ is expanded locally into its Taylor polynomial form, and the contribution of the $k$ -dependent exponential factor is included using numerical algorithms for gamma functions from the Boost.Math library. In order to solve the two resulting square linear systems efficiently, a parallel shared memory direct solver based on sparse LU factorization with pivoting is used, as implemented in the SuperLU library (Demmel, Gilbert & Li Reference Demmel, Gilbert and Li1999; Li et al. Reference Li, Demmel, Gilbert, Grigori, Shao and Yamazaki1999), and the PETSc library (Balay et al. Reference Balay, Adams, Brown, Brune, Buschelman, Eijkhout, Gropp, Kaushik, Knepley, McInnes, Rupp, Smith and Zhang2013) is used for matrix arithmetic. Note that, owing to the dependence of the integral constraint on $k$ , the number of LU factorizations is multiplied by $N_{x}$ . The cost of these factorizations is considerable, and they are tractable only under the condition that $N_{x}$ and $N_{y}$ are not too large.

To cope with the huge scale disparity between the bulk of the channel and the BL, we therefore have to use non-uniform grids in the $y$ direction. During the first phase of the flow evolution, the BL is expected to follow Prandtl’s scaling. The total number of grid points is fixed to $N_{y,1}=385$ . The grid spacing is set to a certain value $\unicode[STIX]{x1D6E5}_{y_{P},1}$ between 0 and

(3.4) $$\begin{eqnarray}\displaystyle L_{1}=166Re^{-1/2}L, & & \displaystyle\end{eqnarray}$$

which corresponds to the BL thickness as can be estimated from the Prandtl calculations, and the remaining points are uniformly spread up to $y=1/2$ , with spacing $\unicode[STIX]{x1D6E5}_{y_{E},1}$ .

Figure 3.

Schematics of discretization grids used by the Navier–Stokes solver for $t\leqslant 54$  (a) and $t>54$  (b). The thickness of the various regions is computed from the Reynolds number as summarized in table 2.

At later times, the Prandtl scaling is expected to break down, and therefore a change of grid is required. For convenience we always perform it at $t=54$ independently of $Re$ . The new grid has $N_{y,2}=449$ points in the $y$ direction. The grid spacing is set to $\unicode[STIX]{x1D6E5}_{y_{K},2}$ between 0 and

(3.5) $$\begin{eqnarray}\displaystyle L_{2}=1210Re^{-1}L, & & \displaystyle\end{eqnarray}$$

then to $\unicode[STIX]{x1D6E5}_{y_{P},2}$ between $L_{2}$ and $L_{1}$ , and the remaining points are spread uniformly up to $y=1/2$ . The values of all deltas are given in table 2, and graphical representations of the two grids are provided in figure 3.

3.3 Prandtl–Euler solver

Following previous work (see e.g. Nguyen van yen, Farge & Schneider Reference Nguyen van yen, Farge and Schneider2009), the Euler equation is approximated by the Navier–Stokes equations with hyper-dissipation, i.e. the dissipation term $\unicode[STIX]{x1D708}\unicode[STIX]{x0394}\unicode[STIX]{x1D714}$ is replaced by $-\unicode[STIX]{x1D708}_{2}(-\unicode[STIX]{x0394})^{2}\unicode[STIX]{x1D714}$ . This approximation is second-order in space (Kato Reference Kato1972), which is sufficient for our purpose here. The boundary conditions are enforced using the classical mirror image technique. Since the vorticity field is antisymmetric with respect to $y=1/2$ , we just need to replace the boundary conditions at $y=0$ and $y=1$ by periodic boundary conditions to effectively impose an exact non-penetration condition. The Navier–Stokes equations are then solved on $\mathbb{T}^{2}$ , taking advantage of the symmetry of the solution, using a fully dealiased sine–cosine pseudo-spectral method corresponding to $N_{x}=N_{y}=4096$ grid points in each direction. A low-storage third-order Runge–Kutta scheme is employed for time discretization, the time step being adjusted dynamically to satisfy the Courant–Friedrichs–Lewy (CFL) condition. The hyperviscosity parameter $\unicode[STIX]{x1D708}_{2}$ was set to $2.146\times 10^{-13}$ , which was found to sufficiently regularize the solution.

To solve the initial value problem for the Prandtl equations (2.16), the spatial domain is first restricted to a finite size $L_{y_{P}}$ in the $y_{P}$ direction, where $L_{y_{P}}$ should be chosen sufficiently large so that the dependence of the solution on its value in the considered time interval is of the same magnitude as other numerical errors. The results presented below were obtained with $L_{y_{P}}=64$ . Spatial discretization is then achieved as for the Navier–Stokes solver, except that the grid in $y$ is regular.

When computing the advection term $v_{P}\unicode[STIX]{x2202}_{y_{P}}\unicode[STIX]{x1D714}_{P}$ , the equations at the edges are obtained by shifting the stencils so that they remain inside the computational grid (no additional condition is required since nonlinear advection vanishes at domain boundaries). For the dissipation term $\unicode[STIX]{x2202}_{y_{P}}^{2}\unicode[STIX]{x1D714}_{P}$ , the integral constraint (2.19) is rewritten as

(3.6) $$\begin{eqnarray}\displaystyle \int _{0}^{\infty }\unicode[STIX]{x2202}_{y_{P}}^{2}\unicode[STIX]{x1D714}_{P}(x,y_{P},t)\,\text{d}y_{P}=\unicode[STIX]{x2202}_{x}p_{E}(x,0,t), & & \displaystyle\end{eqnarray}$$

which is enforced as for the Navier–Stokes solver. Finally, the system is closed by imposing that $\unicode[STIX]{x1D714}_{P}(x,L_{y_{P}},t)=0$ , which is consistent with the fact that the exact solution decays rapidly in $y_{P}$ .

3.4 Interactive boundary layer solver

The interactive solver is similar to the Prandtl solver, the only difference being that the pressure correction given by (2.23) is included at each evaluation of the right-hand side. These modified boundary conditions unfortunately modify the stability region of the time discretization scheme, making it much smaller. We have heuristically derived the constraint

(3.7) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x0394}t<C\frac{\unicode[STIX]{x1D6E5}_{x}^{2}}{\sqrt{\unicode[STIX]{x1D708}}}, & & \displaystyle\end{eqnarray}$$

where $C=1.5$ . For efficiency, we use the non-interactive Prandtl solver up to $t=50$ , and only then do we switch on the interactive term.

3.5 Orr–Sommerfeld solver

To compare the Navier–Stokes solution with the predictions of our linear instability model beyond asymptotics, we have written a simple MATLAB solver for the Orr–Sommerfeld eigenvalue problem. The base velocity profile is taken from the interactive boundary layer computations, and the Orr–Sommerfeld problem is solved independently as desired for each value of $x$ , $k$ and $t$ . The $y_{P}$ variable is again truncated at $L_{y_{P}}=64$ , by using artificial boundary conditions $\unicode[STIX]{x1D714}(L_{y_{P}})=0$ and $\unicode[STIX]{x1D719}^{\prime }(L_{y_{P}})+k\unicode[STIX]{x1D708}^{1/2}\unicode[STIX]{x1D719}(L_{y_{P}})=0$ , which follow from the reconnection with a potential solution at large $y_{P}$ .

A second-order finite difference scheme is used for spatial discretization, written using sparse matrices for efficiency, which leads to a complex, non-symmetric eigenvalue problem. The six eigenvalues with largest imaginary part are solved for using the MATLAB function ‘eigs’, which relies on the implicitly restarted Arnoldi method from ARPACK. As a result, the eigenvalue with the largest imaginary part is readily obtained, and the unstable wavenumber range can thus be detected and estimated.