Hostname: page-component-7c8c6479df-8mjnm Total loading time: 0 Render date: 2024-03-28T10:12:51.956Z Has data issue: false hasContentIssue false

Energy dissipation caused by boundary layer instability at vanishing viscosity

Published online by Cambridge University Press:  26 June 2018

Natacha Nguyen van yen
Affiliation:
Institut für Mathematik, Freie Universität Berlin, Arnimallee 6, 14195 Berlin, Germany
Matthias Waidmann
Affiliation:
Institut für Mathematik, Freie Universität Berlin, Arnimallee 6, 14195 Berlin, Germany
Rupert Klein
Affiliation:
Institut für Mathematik, Freie Universität Berlin, Arnimallee 6, 14195 Berlin, Germany
Marie Farge*
Affiliation:
LMD-CNRS, Ecole Normale Supérieure, 24 rue Lhomond, 75231 Paris CEDEX 5, France
Kai Schneider
Affiliation:
Institut de Mathématiques de Marseille, Aix-Marseille Université and CNRS, Marseille, France
*
Email address for correspondence: marie.farge@ens.fr

Abstract

A qualitative explanation for the scaling of energy dissipation by high-Reynolds-number fluid flows in contact with solid obstacles is proposed in the light of recent mathematical and numerical results. Asymptotic analysis suggests that it is governed by a fast, small-scale Rayleigh–Tollmien–Schlichting instability with an unstable range whose lower and upper bounds scale as $Re^{3/8}$ and $Re^{1/2}$, respectively. By linear superposition, the unstable modes induce a boundary vorticity flux of order $Re^{1}$, a key ingredient in detachment and drag generation according to a theorem of Kato. These predictions are confirmed by numerically solving the Navier–Stokes equations in a two-dimensional periodic channel discretized using compact finite differences in the wall-normal direction, and a spectral scheme in the wall-parallel direction.

Type
JFM Papers
Creative Commons
Creative Common License - CCCreative Common License - BY
This is an Open Access article, distributed under the terms of the Creative Commons Attribution licence (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted re-use, distribution, and reproduction in any medium, provided the original work is properly cited.
Copyright
© 2018 The Authors

1 Introduction

Since the challenge laid down by Euler in 1748 for the Mathematics Prize of the Prussian Academy of Sciences in Berlin, the force exerted onto a solid by a fluid flow has been one of the central unknowns of fluid mechanics. From the vorticity transport equations he had derived, d’Alembert (Reference le Rond d’Alembert1768) deduced his now famous paradox, that this force should vanish, contrary to what experimental results and even everyday observation indicate. A frictional explanation involving the viscosity of the fluid was advanced during the nineteenth century within the frame of the new theories of Navier, Saint-Venant and Stokes, but the actual amplitude of the force remained unaccounted for. Indeed, estimates based on the magnitude of the viscosity $\unicode[STIX]{x1D708}$ of the fluid, or equivalently, after non-dimensionalization, on the inverse of the Reynolds number $Re$ , predict that friction should become negligible when $Re\gg 1$ .

The group working at Göttingen University, notably Prandtl (Reference Prandtl1904) and Blasius (Reference Blasius1908), were first to come up with a way of computing, under some hypotheses, the asymptotic behaviour of the force in the limit $Re\rightarrow \infty$ . Their method relies on the notion that viscous effects are confined to a thin layer along the boundary whose thickness scales like $Re^{-1/2}$ , now called the Prandtl boundary layer (BL), inside which the motion is governed by appropriately rescaled equations, whereas the bulk of the fluid remains inviscid. Momentum transport across this layer gives rise to a net drag force of order $O(Re^{-1/2})$ , which can even be computed explicitly in some academic cases. Although Prandtl’s theory has been very fruitful, it has the drawback of breaking down for sufficiently large $Re$ in practically all relevant flow configurations. Indeed, according to the experiments made in Göttingen itself and elsewhere, the force then acquires the stronger scaling $O(1)$ (Schlichting Reference Schlichting1979).

The precise dynamical mechanism which allows a transition to this regime, starting from a fluid initially at rest, is still unknown today. Prandtl, assuming that the BL approximation becomes invalid beyond a certain separation point along the boundary, had already established that the viscous shear vanishes or, in other words, that the flow in the neighbourhood of the wall reverses direction at this point. After the formal asymptotic expansion achieved by Goldstein (Reference Goldstein1948), which clearly indicates a singularity, a viscous scaling of the parallel coordinate was developed to analyse the flow near the separation point, finally leading to the so-called triple-deck structure (Stewartson Reference Stewartson1974). This steady asymptotic theory is instructive and well developed but remains unsatisfactory, since practically all flows become unsteady above a certain Reynolds number.

The classical understanding about the onset of unsteadiness in shear flows goes back to Rayleigh (Reference Rayleigh1880), who showed in particular that a necessary condition for inviscid instability is that the base velocity profile has an inflection point. Since many BL flows lack such a point, Rayleigh’s result seemed at first to rule out linear instability as a generic mechanism for BL breakdown. This impression was even reinforced by a result of Sommerfeld (Reference Sommerfeld1909) showing that purely linear shear flows ( $u^{\prime \prime }(y)=0$ ) were linearly unconditionally stable even when including the effect of viscosity. But the Göttingen group then realized that including the effects of viscosity with a non-zero second derivative of the base flow could trigger an instability which is absent in the inviscid setting (Prandtl Reference Prandtl1921; Tietjens Reference Tietjens1925). Following some of the ideas developed by Heisenberg (Reference Heisenberg1924) in his thesis, Tollmien (Reference Tollmien1929) then studied asymptotic solutions to the Orr–Sommerfeld equation, and thus achieved an elegant analysis of the instability mechanism. He also produced an approximate marginal curve for what is now known as the Tollmien–Schlichting instability, a widely accepted mechanism for transition to turbulence in BLs. Its range of application is, however, limited to small perturbations around a stationary base flow. The Orr–Sommerfeld equation has also been used to derive rigorously qualitative properties of the solution of the two-dimensional Euler equations – for details we refer to Bedrossian, Masmoudi & Vicol (Reference Bedrossian, Masmoudi and Vicol2016).

Beyond this, one is faced with the full unsteady Navier–Stokes equations (NSEs) at high Reynolds number, and the BL problem becomes so wide ranging that it has been investigated almost independently for the past 60 years by two distinct schools, which we will call the ‘aerodynamical’ and the ‘mathematical’ schools. Aerodynamicists took steady BL theory as their starting point and attempted to generalize the main ideas of Prandtl and Goldstein to the unsteady case. This led in the 1950s to the establishment of the Moore–Rott–Sears criterion (Sears & Telionis Reference Sears and Telionis1975), stating that detachment originates from a point within the BL, not necessarily lying on the wall, where vorticity vanishes and parallel velocity equals that of the exterior flow. Reasoning by analogy with the steady case, Sears & Telionis (Reference Sears and Telionis1975) conjectured that separation coincides with the appearance of a singularity in the solution at some finite time $t^{\ast }$ . The approach to such a singularity was confirmed using numerical experiments by Van Dommelen & Shen (Reference Van Dommelen and Shen1980) for the impulsively started cylinder, and then analysed in detail by power-series expansions in Lagrangian variables by Van Dommelen & Cowley (Reference Van Dommelen and Cowley1990). We refer the reader to Van Dommelen’s (Reference Van Dommelen1981) PhD thesis for a detailed pedagogical account of these findings. Later work largely supports his results (Ingham Reference Ingham1984; Weinan & Engquist Reference E and Engquist1997; Gargano, Sammartino & Sciacca Reference Gargano, Sammartino and Sciacca2009).

The next challenge was to understand what happens to the actual NSE flow as the corresponding Prandtl solution becomes singular. Elliott, Smith & Cowley (Reference Elliott, Smith and Cowley1983) obtained the estimate $(t^{\ast }-t)=O(Re^{-2/11})$ for the time $t$ at which the BL assumptions first break down. To describe the solution at later times, some hope came from the interacting boundary layer (IBL) method, which relieves the Goldstein singularity in steady BLs, but it was quickly shown (Smith Reference Smith1988; Peridier, Smith & Walker Reference Peridier, Smith and Walker1991a ,Reference Peridier, Smith and Walker b ) to lead to a finite-time singularity when applied to unsteady problems.

Over the same period of time, the mathematical school focused on totally different issues concerning the initial value problem for the unsteady Prandtl equations, on the one hand, and the vanishing viscosity limit for solutions of the NSEs on the other. Local well-posedness for the Prandtl equations was proved long ago by Oleinik (Reference Oleinik1966) in cases where detachment is not expected, i.e. for monotonic initial data and favourable pressure gradient. Sammartino & Caflisch (Reference Sammartino and Caflisch1998) showed local well-posedness without the monotonicity conditions, but this time under a very harsh regularity condition, that the amplitude of the parallel Fourier coefficients of all quantities decrease exponentially with wavenumber. These conditions were recently improved by Gérard-Varet & Masmoudi (Reference Gérard-Varet and Masmoudi2015), but the required decay of Fourier coefficients is still faster than any power of $k$ . Although no rigorous construction of such a solution exists to our knowledge, it is generally believed that there exist exact solutions of Prandtl equations which blow up in finite time. Maekawa (Reference Maekawa2014) proved the convergence of the Euler and Prandtl solution versus the Navier–Stokes (NS) solution in the $L^{\infty }$ norm with order $\sqrt{\unicode[STIX]{x1D708}}$ for an initial condition where the vorticity is compactly supported at finite distance from the wall.

Kato (Reference Kato1984) made a decisive contribution, by linking the vanishing viscosity limit problem to the behaviour of energy dissipation at the boundary. Kato’s theorem implies, in the particular case of a flow in a smooth two-dimensional domain $\unicode[STIX]{x1D6FA}$ with smooth initial data and without forcing, the equivalence between the following assertions:

  1. (i) The NS flow converges to the Euler flow uniformly in time in the energy norm.

  2. (ii) The energy dissipation associated with the NS flow, integrated over a strip of thickness proportional to $Re^{-1}$ around the solid, which we will call the Kato layer, tends to zero as $Re\rightarrow \infty$ .

Since convergence to the Euler flow excludes detachment, one of the essential messages carried by this theorem is that the flow has to develop dissipative activity at a scale at least as fine as $Re^{-1}$ for detachment to be possible. Later refinements of Kato’s work linked breakdown to scalings $O(Re^{\unicode[STIX]{x1D6FC}})$ with $\unicode[STIX]{x1D6FC}\geqslant 3/4$ of the wall pressure gradient (Temam & Wang Reference Temam and Wang1997), and $O(Re^{\unicode[STIX]{x1D6FD}})$ with $\unicode[STIX]{x1D6FD}\geqslant 1/2$ of the $L^{2}$ norm of velocity in the Kato layer (Kelliher Reference Kelliher2007). For a detailed review about mathematical theorems related to turbulence, we refer to Bardos & Titi (Reference Bardos and Titi2013).

Not only is the gap between the aerodynamical and the mathematical schools of thought quite impressive when one realizes that they are really concerned with the same problem, but after close consideration, it even appears that they contradict each other on an essential point, which we now attempt to clarify. As shown by Van Dommelen & Cowley (Reference Van Dommelen and Cowley1990), the finite-time singularity in the unsteady Prandtl equations, which is characterized by a blow-up of parallel vorticity gradients, does correspond to a ‘detachment’ process, in the sense that there exist fluid particles that are accelerated infinitely rapidly away from the wall. In the following, we shall call this process ‘eruption’. Since its discovery, it seems to have been at least tacitly assumed to underlie the initial stage of the (a priori different) detachment process actually experienced by the NS solution. According to this scenario, singularity would be avoided in the NS case thanks to a large-scale process not taken into account in the Prandtl approximation, namely the normal pressure gradient. The Kato criterion, on the other hand, tells us something entirely different, which is that, for detachment to happen, scales as fine as $Re^{-1}$ , which are not even accounted for in the Prandtl solution, need to come into play. This suggests that eruption never really enters the scene in the NS flow, being indeed short-circuited by a faster mechanism at finer scales.

In the last decade, instabilities at high parallel wavenumber came up as a possible explanation for these finer scales. On the mathematical side, Grenier (Reference Grenier2000) proved that Prandtl’s asymptotic expansion is invalid for some types of smooth perturbed shear flows, due to instabilities at high parallel wavenumbers. Then Gérard-Varet & Dormy (Reference Gérard-Varet and Dormy2010) showed, again for smooth perturbed shear flows, that the Prandtl equations could be linearly ill-posed in any Sobolev space (i.e. assuming spectra which decay like powers of $k$ ), although they are locally well-posed in the analytical framework, as we have seen above. In the last decade, Grenier, Guo & Nguyen (Reference Grenier, Guo and Nguyen2015, Reference Grenier, Guo and Nguyen2016) have continued working on the ambitious programme aiming to achieve a rigorous mathematical description of instabilities in generic BL flows. Once again, their proofs show that the ill-posedness is due to modes with large wavenumber in the direction parallel to the wall. Grenier & Nguyen (Reference Grenier and Nguyen2018) even studied higher-order terms in the asymptotic expansion of the inviscid limit but were faced with the appearance of further instabilities which have not been tamed up to now.

Several fluid dynamicists have also advanced the idea that instability-type mechanisms may play an important role in the process of unsteady detachment. Cassel (Reference Cassel2000) directly compared numerical solutions of the NSEs and of the corresponding Prandtl equations in an attempt to verify the correctness of the asymptotic expansion. Although he considered a vortex-induced BL instead of the impulsively started cylinder, his Prandtl solution behaves qualitatively like the one in Van Dommelen & Shen (Reference Van Dommelen and Shen1980), developing strong parallel velocity gradients in a process which seems to lead to a finite-time singularity associated with a large normal displacement of fluid particles concentrated around a single parallel location. But interestingly, around the same time, the NS solution adopts a quite different behaviour, characterized by the appearance of strong oscillations in the wall-parallel pressure gradient, which he was not able to explain. Brinckman & Walker (Reference Brinckman and Walker2001) also saw oscillations, and claimed that they were due to a Rayleigh instability of the shear velocity profile, a hypothesis which was developed further in the review paper by Cowley (Reference Cowley, Aref and Phillips2002). Although the numerics underlying these findings were later invalidated by Obabko & Cassel (Reference Obabko and Cassel2002) due to their insufficient grid resolution, the existence of an instability mechanism has continued to be a hot topic in later papers on unsteady detachment (Bowles, Davies & Smith Reference Bowles, Davies and Smith2003; Bowles Reference Bowles2006; Cassel & Obabko Reference Cassel and Obabko2010; Gargano, Sammartino & Sciacca Reference Gargano, Sammartino and Sciacca2011; Gargano et al. Reference Gargano, Sammartino, Sciacca and Cassel2014). In fact, it seems to be the only surviving conjecture at this time to explain unsteady detachment. The nature and quantitative properties of the instability remain to be elucidated.

Imposing no-slip boundary conditions on high-Reynolds-number flows, even in two dimensions, is a tough numerical problem and one should be especially careful with the numerics given that the problem is theoretically not well understood yet. In our previous work (Nguyen van yen, Farge & Schneider Reference Nguyen van yen, Farge and Schneider2011), we computed a series of dipole–wall collisions, a well-studied academic flow introduced by Orlandi (Reference Orlandi1990). Our goal was to derive the scaling of energy dissipation when the Reynolds number increases, for fixed initial data and geometry. According to the Kato criterion, this is an important element to understand detachment. We chose to work with a volume penalization scheme, which has the advantages of being efficient, easy to implement and, most importantly, providing good control on numerical dissipation. However, the no-slip condition is only approximated, and the higher the Reynolds number, the more costly it is to enforce satisfactorily. In fact, post-processing numerically calculated flows revealed that they effectively experienced Navier boundary conditions with a slip length proportional to $Re^{-1}$ – for a more detailed analysis of the scheme, see Nguyen van yen, Kolomenskiy & Schneider (Reference Nguyen van yen, Kolomenskiy and Schneider2014). In this setting, we did find indications that energy dissipation converges to a finite value when $Re\rightarrow \infty$ .

Sutherland, Macaskill & Dritschel (Reference Sutherland, Macaskill and Dritschel2013) reconsidered the computations of Nguyen van yen et al. (Reference Nguyen van yen, Farge and Schneider2011) and studied the effects of a finite slip length inversely proportional to the Reynolds number. They confirmed our findings using compact finite differences on a Chebyshev grid in the wall-normal direction and a Fourier method in the wall-parallel direction. Therewith exact no-slip boundary conditions and Navier boundary conditions with a given slip length (independently of $Re$ ) can be imposed. In the no-slip case and also for a fixed slip length, they concluded that energy dissipation vanishes and that there is no indication of the persistence of energy-dissipating structures in the vanishing viscosity limit. This is quite unexpected since, due to the spectral properties of the Stokes operator, the no-slip boundary condition is stiffer that any Navier boundary condition with a non-zero slip length, and should thus generate larger gradients, or in other words more dissipation. This makes the claim of Sutherland et al. (Reference Sutherland, Macaskill and Dritschel2013), that there is dissipation at vanishing viscosity with Navier conditions for a fixed slip length, but not with no-slip conditions, quite counterintuitive. In fact, looking more closely at their results, it appears that the central claim is based on a single computation (see figures 17 and 18 of their paper), for which a convergence test is not provided, which in our opinion leaves the matter unsettled.

Clercx & van Heijst (Reference Clercx and van Heijst2002) already studied the role of no-slip boundaries in two-dimensional flows considering the dipole–wall interaction using a Chebyshev spectral method. They observed enstrophy scalings $Z\propto Re^{0.8}$ for $Re\leqslant 20\,000$ and $Z\propto Re^{0.5}$ for $Re\geqslant 20\,000$ , where $Re$ is based on the velocity and the size of the dipole. Then Keetels et al. (Reference Keetels, Kramer, Clercx and van Heijst2011) proposed an oscillating plate model as an analogous simplified boundary layer to predict these scalings, which are in reasonable agreement with the numerical simulations of Clercx & van Heijst (Reference Clercx and van Heijst2002). These scaling observations indeed imply that dissipation would vanish in the limit of infinite Reynolds number. However, dissipation effects remain highly important for high-Reynolds-number two-dimensional flows in wall-bounded domains, as reviewed recently by Clercx & van Heijst (Reference Clercx and van Heijst2017).

Following the findings of Sutherland et al. (Reference Sutherland, Macaskill and Dritschel2013), we decided to revisit the issue once again, but this time using a numerical scheme which has both high precision and accurate no-slip boundary conditions. For this, we have turned to compact finite differences with an ad hoc irregular grid in the wall-normal direction, and Fourier coefficients in the wall-parallel direction. Combining ideas developed in the last decade by many authors, we propose a heuristic scenario for detachment, based on an instability mechanism of the Tollmien–Schlichting type, which also explains the new vorticity scaling $O(Re^{1})$ and the occurrence of dissipation. We then check that all these processes are actually occurring in our numerical solution of the dipole–wall initial value problem. For this, adopting a methodology similar to the one employed by Cassel (Reference Cassel2000), we directly compare the NS solution and the corresponding Prandtl approximation which we also compute.

In the next section (§ 2), we introduce the flow configuration and the corresponding NS and Prandtl models. Although these models are classical, we present specific reformulations which were chosen in order to facilitate both numerical efficiency and theoretical interpretation. Then we use the model to predict the appearance of an instability and understand its characteristics. In § 3, we introduce the model discretizations which we have used for our numerical computations. In § 4 we describe the numerical results. Finally, in § 5 we analyse the numerical results in the light of our preceding theoretical developments and we draw the necessary conclusions.

2 Models

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}$$
which in turn satisfies the Poisson equation
(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}$$
where $p_{E}$ is the pressure field calculated from $\unicode[STIX]{x1D714}_{E}$ . It is instructive to rederive the classical Neumann condition (2.16d ) as follows. First, by replacing $\unicode[STIX]{x1D714}$ according to (2.15) in (2.12), one obtains
(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 Solvers

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}$$
where the coefficients are calculated by matching the Taylor expansions of both sides up to fifth order.

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.

4 Results

4.1 Before detachment

The behaviour of the various solutions well before the Prandtl singularity time is well understood and we present it only for the sake of completeness. After a rapid relaxation phase, the initial vorticity distribution splits into two counter-propagating dipoles, each of which shoots towards one of the channel walls. At this point, the Navier–Stokes vorticity field in the bulk flow remains similar to the Euler vorticity field, as shown by comparing their contour lines at $t=30$ in figure 4(a). The Navier–Stokes flow in the Prandtl boundary layer units (figure 4 c) is smooth, and well approximated by the corresponding solution of the Prandtl equations, shown in blue. As the dipole approaches the wall, the pressure gradient increases, causing inward diffusion of vorticity as well as increased vorticity gradients within the boundary layer. At $t=50$ , we still observe qualitative similarity between the Navier–Stokes flow at high Reynolds number, on the one hand, and the Euler flow in the bulk with the Prandtl flow in the boundary layer, on the other (figure 4 b,d). However, as expected, the discrepancy between Prandtl and Navier–Stokes flows has increased.

Figure 4. Comparison between Navier–Stokes at $Re=123\,075$ (black) and Euler/Prandtl (blue). (a,b) Contour lines of vorticity field at $t=30$  (a) and $t=50$  (b), for $\unicode[STIX]{x1D714}=0.1,0.2,\ldots ,1$ . (c,d) Contour lines of corresponding boundary layer vorticity fields in Prandtl variables, for $\unicode[STIX]{x1D714}=-|\unicode[STIX]{x1D714}|_{max}/6,-2|\unicode[STIX]{x1D714}|_{max}/6,\ldots ,-5|\unicode[STIX]{x1D714}|_{max}/6$ , where $|\unicode[STIX]{x1D714}|_{max}=3.75\times 10^{-4}$ and $|\unicode[STIX]{x1D714}|_{max}=8.70\times 10^{-3}$ , respectively.

Figure 5. Contour lines of vorticity field in Prandtl variables at $t=54$  (a) and $t=55.3$  (b), for $\unicode[STIX]{x1D714}=5|\unicode[STIX]{x1D714}|_{max}/6,\ldots ,|\unicode[STIX]{x1D714}|_{max}/6,0,-|\unicode[STIX]{x1D714}|_{max}/6,-2|\unicode[STIX]{x1D714}|_{max}/6,\ldots ,-5|\unicode[STIX]{x1D714}|_{max}/6$ , where $|\unicode[STIX]{x1D714}|_{max}=1.25\times 10^{-2}$ and $|\unicode[STIX]{x1D714}|_{max}=1.30\times 10^{-2}$ , respectively. Black contour lines correspond to Navier–Stokes at $Re=123\,075$ , while blue lines correspond to Euler/Prandtl. The $\unicode[STIX]{x1D714}=0$ contour line is shown in bold.

At $t=54$ (figure 5 a) a new important feature of the flow is that a region of opposite-sign vorticity has appeared within the boundary layer, indicating the build-up of a recirculation bubble along the wall. This effect is well captured by the Prandtl flow, which overall continues to approach the Navier–Stokes solution pretty well, although the discrepancy has again notably increased.

4.2 Prandtl blow-up

The first signs of a qualitatively different behaviour become visible shortly thereafter, as shown for example at $t=55.3$ in figure 5(b). The contour lines of the Prandtl vorticity have become very concentrated around $x^{\ast }=0.556$ , indicating the formation of a finite-time singularity with precisely the qualitative features predicted by Van Dommelen & Cowley (Reference Van Dommelen and Cowley1990), in particular a blow-up of the wall-normal velocity associated with an abrupt acceleration of fluid particles away from the wall.

As the Prandtl solution approaches its singularity time $t^{\ast }\simeq 55.6$ , parallel vorticity gradients increase rapidly, and soon the cut-off parallel wavenumber of the numerical scheme becomes insufficient to resolve it.

4.3 Large-scale interaction and instability

According to Elliott et al. (Reference Elliott, Smith and Cowley1983), the Navier–Stokes solution departs from the Prandtl behaviour when the potential flow perturbation due to the presence of the boundary layer starts to perturb the wall pressure gradient.

Figure 6. Comparison between Navier–Stokes (black), Prandtl–Euler (blue) and interactive boundary layer (red) models at $t=50$ (a) and $t=54$ (b). Contour lines of vorticity (a,b) and vorticity along the boundary (c,d). The $\unicode[STIX]{x1D714}=0$ contour line is shown in bold.

By comparing the Navier–Stokes, Prandtl–Euler and interactive boundary layer models at different Reynolds numbers for $t=50$ , we observe good agreement between all models (figure 6 a). For $t=54$ (figure 6 b) it can be noted that the IBL solution has indeed slightly departed away from the Prandtl solution, but to a point which falls way short of capturing the full behaviour of the NS solution. This effect, which corresponds in principle to the large-scale interaction also described by Gargano et al. (Reference Gargano, Sammartino, Sciacca and Cassel2014), seems to play only a secondary role in our setting. More importantly, we observe the growth of an elbow feature in the NS solution, indicating the start of the growth of a packet of higher- $k$ modes concentrated around $x=0.6$ .

Figure 7. Vorticity along the boundary for Navier–Stokes (solid line) for varying Reynolds number (from $Re=7692$ to $123\,075$ ) and Prandtl (dashed line).

Figure 8. Fourier spectra of boundary vorticity for Navier–Stokes (solid line) at different Reynolds numbers (from $Re=7692$ to $123\,075$ ) and Prandtl (dashed line). The spectra have been smoothed to eliminate fast oscillations due to the odd symmetry of the function (see text).

When considering the NSE solution for varying $Re$ at the instants $t=54$ and $t=55.3$ in figure 7, we observe that for increasing $Re$ the elbow structure looks more and more like a wavepacket confined to a well-defined interval on the wall. To understand better the onset of these oscillations, it is tempting to consider one-dimensional Fourier transforms of those wall vorticity traces. Unfortunately, the odd symmetry of the function around the dipole axis gives rise to fast oscillations in the Fourier coefficients which impair the readability of the spectra. To get rid of this effect, the spectra are averaged out using the low-pass filter $\exp (-(5.2|k|/N)x^{16})$ . The results are shown in figure 8.

The Prandtl solution (dashed curves) develops a $k^{-3/2}$ power-law profile at high $k$ , consistent with the build-up of a jump singularity in $\unicode[STIX]{x1D714}$ along the wall. In contrast, the NSE solution spectra decay exponentially at high $k$ and develop a distinctive bump in the wavenumber range. Both the width and $k$ location of this bump increase with Reynolds number and with time. Interestingly, for the largest $Re$ considered, a relatively good separation of scales can be observed at $t=55.3$ between the low- $k$ features and the high- $k$ wavepacket, the transition occurring somewhere between $k=20$ and $k=30$ . This confirms a posteriori the validity of the slowly varying flow approximation used in deriving the asymptotic stability results of § 2.4. Moreover, all solutions have exponentially decaying spectra at sufficiently large values of $k$ , consistent with their analytic regularity being well resolved in the current numerical setting.

Figure 9. Stability boundary subject to the frozen flow approximation, shown as a function of $x$ at $t=54$ and for $Re=123\,075$ . The upper and lower marginal wavenumbers estimated from our asymptotic analysis (2.51) and (2.54) are shown in dashed lines.

We would now like to compare the characteristics of these spectra with the predictions of our analysis based on the Orr–Sommerfeld equations. The numerical instability range obtained by direct eigenmode analysis (figure 9, bold line) extends over the interval $x\in [0.555,0.65]$ on the boundary, which is in good qualitative agreement with the spatial extent of the oscillations seen in figure 6(b). In $k$ -space, the Orr–Sommerfeld computations predict that the instability should start around $k=30$ for the high $Re$ considered, which is in very good agreement as well with the wavenumber at which the corresponding spectrum in figure 8(a) starts to exceed the reference Prandtl solution (shown in dashed lines). This effect becomes more pronounced at $t=55.3,$ shown in figure 8(b). Another important point consistent with our scenario is that the stable modes $10<k<30$ indeed appear damped in the NSE solutions compared to the Prandtl solution, a phenomenon which would be very hard to explain using a singularity-type scenario.

Concerning the theoretical prediction for the lower end of the instability range, qualitative agreement is restricted to a narrow region around $x=0.6$ , whereas the wavenumber is very underestimated as soon as a point where $u_{P}^{\prime }(0)=0$ is approached. Nevertheless, the overall instability region is qualitatively well captured by the criterion $u_{P}^{\prime }(0)<0$ .

4.4 Detachment and production of dissipative structures

Figure 10. Contour lines of vorticity field in Prandtl variables for Navier–Stokes at $t=56$  (a) and $t=56.9$  (b). The $\unicode[STIX]{x1D714}=0$ contour line is shown in bold.

Figure 11. Vorticity along the boundary for various Reynolds numbers (from $Re=7692$ to $123\,075$ ) in Prandtl variables at $t=56$ (a) and $t=56.9$ (b).

The instability process which we have seen at play above introduces a new vorticity scaling, $Re^{1}$ , very close to the wall. This new scaling is difficult to notice at first, because it is hidden behind the pre-existing large negative vorticity of the boundary layer. A simple trick to observe it more easily is to consider only the maximum vorticity of positive sign (figure 12 a). This quantity scales like $Re^{1/2}$ at $t=53$ , and at $t=56.9$ it has clearly transited to the stronger $Re^{1}$ scaling. Accordingly, the enstrophy scaling has become dissipative at $t=56.9$ , thus indicating the production of a dissipative structure as predicted by the Kato criterion. Shortly thereafter, several further extrema with alternating signs successively appear for sufficiently high Reynolds number, corresponding to increasingly fine parallel scales in the $x$ direction, as illustrated in figures 10 and 11, and previously observed in figure 12 of Kramer, Clercx & van Heijst (Reference Kramer, Clercx and van Heijst2007). Note that in table V of Kramer et al. (Reference Kramer, Clercx and van Heijst2007) the vorticity maxima have been reported for $625\leqslant Re\leqslant 20\,000$ . Plotting their values as a function of Reynolds number yields scalings from $Re^{1}$ to $Re^{0.65}$ .

Figure 12. Reynolds-number scalings of maximum vorticity (a) and enstrophy (b) for Prandtl and Navier–Stokes before detachment ( $t=53$ ), and for Navier–Stokes after detachment ( $t=56.9$ ).

4.5 Later evolution

Figure 13. Contour lines of vorticity field at $t=80$ for Navier–Stokes at $Re=15\,384$ and Euler/Prandtl solutions.

At much later times, the Euler and Navier–Stokes solutions have become entirely different (figure 13). In the Euler case, the vortices glide along the wall, having paired with their mirror image, and no new vorticity has been generated. Energy and enstrophy are both conserved. In the Navier–Stokes case, the detachment process has led to the formation of two new vortices, of much larger amplitude than those of the incoming dipole. The activity in the boundary layer remains intense, leading to the ejection of smaller structures.

4.6 Convergence checks

An essential point concerns the control of the discretization error. Following common practice in numerical fluid dynamics, we have taken care to use quite pessimistic scalings to design the wall-parallel and wall-normal grids in order to resolve the necessary range of scales, and as a result we have not observed spurious grid-scale oscillations which would suffice to indicate under-resolution.

Figure 14. Comparison between Navier–Stokes solutions at the same Reynolds number ( $Re=123\,075$ ) with two numerical resolutions, in order to check for convergence: (a) $t=56$ , (b) $t=56.9$ .

To go one step further, we now consider the contour lines of vorticity during detachment, at $t=56$ and $t=56.9$ . As shown in figure 14, the computation used in our analysis and the one obtained with twice the grid spacing agree well.

5 Discussion

Several features of our numerical solutions had not been observed in previous work. The most striking one is the appearance of the scaling $Re^{1}$ for the vorticity maximum, which takes precedence, at the singularity time, over the weaker Prandtl scaling $Re^{1/2}$ . Even more strikingly, as seen on the graphs of wall vorticity, this new extremum does not even appear at the same location as the Prandtl singularity. This result contradicts sharply the picture of boundary layer detachment as it was described in earlier work, as an essentially local process coinciding with the singularity in the Prandtl equations. Thanks to the vorticity formulation we have favoured, the origin of the non-locality can be traced back to the integral constraints (2.12) on the vorticity field, which are themselves consequences of the no-slip boundary conditions combined with incompressibility. If higher and higher $k$ modes are excited, as occurs in particular due to the Prandtl singularity formation, the reaction of the flow dictated by (2.1) has no reason to be localized in the $x$ direction.

A key observation is that, in regions with reversed flow near the wall, the width of the unstable wavenumber range scales like $Re^{1/2}$ , while the amplitude of vorticity continues to scale as $Re^{1/2}$ due to the presence of a Prandtl boundary layer. Therefore, as soon as the build-up of the Prandtl singularity sufficiently excites those wavenumbers, their superposition induces a $Re^{1}$ scaling for the amplitude of $\unicode[STIX]{x1D714}$ . In the linear phase, the thickness of the corresponding new wall-normal sublayer scales like $Re^{-2/3}$ , but as soon as the instability becomes nonlinear, vorticity transport induces excitation of scales as fine as $Re^{-1}$ , leading to dissipation. According to this scenario, the process of detachment is thus intricately linked to the occurrence of dissipation.

Another open question concerns the description of the flow after detachment. If it is confirmed, the scenario we are proposing indicates that the process of detachment and vorticity production by no-slip walls could be modelled by detecting Prandtl singularities and, when they are about to occur, by introducing nonlinear Rayleigh–Tollmien–Schlichting waves, followed by roll-up and the injection of a dissipative structure into the bulk flow. However, an essential point to keep in mind is that the phase of these waves is very sensitive to Reynolds number, which means that there is little hope of a fully deterministic Reynolds-number-independent description. This could have important consequences for the modelling of wall-bounded turbulent flows.

The existence of vortical structures in turbulent boundary layers is well established (Robinson Reference Robinson1991). The local conditions in such flows are therefore not as different from those we have studied as one might first expect. According to the logarithmic law of the wall

(5.1) $$\begin{eqnarray}\displaystyle \langle U(y)\rangle \simeq \frac{U_{\unicode[STIX]{x1D70F}}}{K_{karman}}\log \left(\frac{yU_{\unicode[STIX]{x1D70F}}}{\unicode[STIX]{x1D708}}\right), & & \displaystyle\end{eqnarray}$$

where

(5.2) $$\begin{eqnarray}\displaystyle U_{\unicode[STIX]{x1D70F}}=\sqrt{\unicode[STIX]{x1D708}\left\langle \left.\frac{\text{d}U}{\text{d}y}\right|_{y=0}\right\rangle } & & \displaystyle\end{eqnarray}$$

is the so-called friction velocity. This behaviour is confirmed by the most recent experiments, with subtle corrections (see e.g. Marusic et al. Reference Marusic, McKeon, Monkewitz, Nagib, Smits and Sreenivasan2010). An important consequence is that the bulk velocity and $U_{\unicode[STIX]{x1D70F}}$ have the same scaling with $Re$ up to a logarithmic factor. Then, from (5.2), one can immediately deduce that $\text{d}\langle u\rangle /\text{d}y|_{y=0}$ scales like $Re^{1}$ up to a logarithmic factor, which can be seen as the statistical signature of the existence of a boundary layer of thickness $Re^{-1}$ in the neighbourhood of the wall. Hence we see that the log law, as an experimental result, is consistent in some sense with the existence of a Kato layer, as we have established in our two-dimensional computations in a much more restricted setting. This connection can be made, as we just did, in a purely phenomenological way without invoking the Kolmogorov scale and the notion of cascade. In fact, the essential point is that $U_{\unicode[STIX]{x1D70F}}$ scales with the bulk velocity, and this scaling was introduced by von Kármán (Reference von Kármán1921) precisely to account for the behaviour of the drag coefficient at high Reynolds number, which was recognized as the essential issue at this time. From this discussion it appears that our results may help in investigating rigorous foundations to the phenomenological Kármán theory.

Acknowledgements

The authors would like to thank C. Bardos, D. Gérard-Varet, H. N. Lopes, M. L. Filho and especially S. Cowley for important discussions. N.N.V.Y. thanks the Humboldt Foundation for supporting his research by a post-doctoral fellowship. M.F. and K.S. acknowledge support by the French Research Federation for Fusion Studies within the framework of the European Fusion Development Agreement (EFDA).

Supplementary movie

A supplementary movie is available at https://doi.org/10.1017/jfm.2018.396.

Appendix A. Estimation of the integral $I_{c}$

In this appendix, we establish the estimate (2.47) for the integral $I_{c}$ defined by (2.43). For simplicity we drop the index $P$ , and denote generically by $K_{R}$ real numbers not depending on $c$ . Our guess is that the dominant behaviour of $I_{c}$ when $c\rightarrow 0$ is controlled by the behaviour of $u$ around its zeros on $[0,+\infty [$ . Hence, let $(y_{i})_{1\leqslant i\leqslant m}$ denote the locations of these zeros. We assume that $u^{\prime }(y_{i})\neq 0$ for all $i$ , and that $u$ has an analytic continuation in a complex neighbourhood of the real axis, which ensures that there exist complex neighbourhoods $Y_{i}$ of $y_{i}$ such that $u$ is a local holomorphism $u:Y_{i}\rightarrow u(Y_{i})$ , and $u(Y_{i})$ are neighbourhoods of $0$ . Assuming that $c$ is sufficiently close to zero so that it lies in the intersection of all the $u(Y_{i})$ and it is smaller than the infimum of $u$ outside of the $Y_{i}$ , the equation $c=u(z)$ has exactly $m$ solutions, which we denote $z_{i}(c)$ . Now letting

(A 1) $$\begin{eqnarray}\displaystyle f_{i}(y,c)=\frac{1}{(u(y)-c)^{2}}-\frac{1}{(u_{\infty }-c)^{2}}-\frac{1}{u^{\prime }(z_{i}(c))^{2}(y-z_{i}(c))^{2}}+\frac{u^{\prime \prime }(z_{i}(c))}{u^{\prime }(z_{i}(c))^{3}(y-z_{i}(c))}, & & \displaystyle \nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$

it follows from a straightforward Taylor expansion of $u$ around $z_{i}(c)$ that $f_{i}(y,c)$ is bounded on $Y_{i}$ by a constant independent of $c$ , and therefore equi-integrable on $Y_{i}$ . Therefore, if we define $\unicode[STIX]{x1D709}_{0}=0$ , $\unicode[STIX]{x1D709}_{i}=(y_{i+1}+y_{i})/2$ for $i<m$ and $\unicode[STIX]{x1D709}_{m}=2y_{m}+1$ , $f_{i}(y,c)$ is equi-integrable on $[\unicode[STIX]{x1D709}_{i},\unicode[STIX]{x1D709}_{i+1}]$ , and by the Vitali convergence theorem, the pointwise limit

(A 2) $$\begin{eqnarray}\displaystyle f_{i}(y,0)=\frac{1}{u(y)^{2}}-\frac{1}{u^{\prime }(y_{i})^{2}(y-y_{i})^{2}}+\frac{u^{\prime \prime }(y_{i})}{u^{\prime }(y_{i})^{3}(y-y_{i})}-\frac{1}{u_{\infty }^{2}} & & \displaystyle\end{eqnarray}$$

obtained when $c\rightarrow 0$ is integrable, its integral being the limit of the integrals of $f_{i}(y,c)$ when $c\rightarrow 0$ , i.e.

(A 3) $$\begin{eqnarray}\displaystyle \int _{\unicode[STIX]{x1D709}_{i}}^{\unicode[STIX]{x1D709}_{i+1}}f_{i}(y,c)\,\text{d}y\xrightarrow[{c\rightarrow 0}]{}\int _{\unicode[STIX]{x1D709}_{i}}^{\unicode[STIX]{x1D709}_{i+1}}f_{i}(y,0)\,\text{d}y. & & \displaystyle\end{eqnarray}$$

Now letting

(A 4a ) $$\begin{eqnarray}\displaystyle & \displaystyle J_{i,c}=\int _{\unicode[STIX]{x1D709}_{i}}^{\unicode[STIX]{x1D709}_{i+1}}\frac{\text{d}y}{u^{\prime }(z_{i}(c))^{2}(y-z_{i}(c))^{2}}, & \displaystyle\end{eqnarray}$$
(A 4b ) $$\begin{eqnarray}\displaystyle & \displaystyle K_{i,c}=-\int _{\unicode[STIX]{x1D709}_{i}}^{\unicode[STIX]{x1D709}_{i+1}}\frac{u^{\prime \prime }(z_{i}(c))\,\text{d}y}{u^{\prime }(z_{i}(c))^{3}(y-z_{i}(c))}, & \displaystyle\end{eqnarray}$$
by construction
(A 5) $$\begin{eqnarray}I_{c}=-c\mathop{\sum }_{i=0}^{m-1}\left(\int _{\unicode[STIX]{x1D709}_{i}}^{\unicode[STIX]{x1D709}_{i+1}}f_{i}(y,c)\,\text{d}y+J_{i,c}+K_{i,c}\right)-c\int _{\unicode[STIX]{x1D709}_{m}}^{\infty }\left(\frac{1}{(u(y)-c)^{2}}-\frac{1}{(u_{\infty }-c)^{2}}\right)\end{eqnarray}$$

or, by (A 3) and the fact that $u$ does not vanish on $[\unicode[STIX]{x1D709}_{m},+\infty ]$ ,

(A 6a ) $$\begin{eqnarray}\displaystyle I_{c} & \underset{c\rightarrow 0}{=} & \displaystyle \displaystyle -c\mathop{\sum }_{i=0}^{m-1}\left(\int _{\unicode[STIX]{x1D709}_{i}}^{\unicode[STIX]{x1D709}_{i+1}}f_{i}(y,0)\,\text{d}y+J_{i,c}+K_{i,c}\right)-c\int _{\unicode[STIX]{x1D709}_{m}}^{\infty }\left(\frac{1}{u(y)^{2}}-\frac{1}{u_{\infty }^{2}}\right)+o(c)\qquad\end{eqnarray}$$
(A 6b ) $$\begin{eqnarray}\displaystyle & \underset{c\rightarrow 0}{=} & \displaystyle \displaystyle -c\mathop{\sum }_{i=0}^{m-1}(J_{i,c}+K_{i,c})+cK_{R}+o(c)\end{eqnarray}$$
so that the dominant behaviour of $I_{c}$ is controlled by $J_{i,c}$ and $K_{i,c}$ , which we now proceed to estimate. Since
(A 7) $$\begin{eqnarray}\displaystyle J_{i,c}=-\frac{1}{u^{\prime }(z_{i}(c))^{2}}\left(\frac{1}{\unicode[STIX]{x1D709}_{i+1}-z_{i}(c)}-\frac{1}{\unicode[STIX]{x1D709}_{i}-z_{i}(c)}\right), & & \displaystyle\end{eqnarray}$$

$J_{i,c}$ converges to a real constant $J_{i}$ if $i\neq 0$ , but diverges for $i=0$ . By Taylor-expanding $u$ around $0$ , we get that

(A 8) $$\begin{eqnarray}\displaystyle z_{0}(c)\underset{c\rightarrow 0}{=}\frac{c}{u^{\prime }(0)}+o(c), & & \displaystyle\end{eqnarray}$$

and therefore

(A 9) $$\begin{eqnarray}\displaystyle J_{0,c}\underset{c\rightarrow 0}{=}-\frac{1}{u^{\prime }(0)c}+K_{R}+o(1). & & \displaystyle\end{eqnarray}$$

Now

(A 10) $$\begin{eqnarray}\displaystyle K_{i,c}=-\frac{u^{\prime \prime }(z_{i}(c))}{u^{\prime }(z_{i}(c))^{3}}[\ln (\unicode[STIX]{x1D709}_{i+1}-z_{i}(c))-\ln (\unicode[STIX]{x1D709}_{i}-z_{i}(c))], & & \displaystyle\end{eqnarray}$$

where for the complex logarithm we have legitimately taken the principal branch, since the integration path does not cross the negative real axis. As before we first assume that $i\neq 0$ , in which case

(A 11a ) $$\begin{eqnarray}\displaystyle K_{i,c} & \underset{c\rightarrow 0}{=} & \displaystyle -\frac{u^{\prime \prime }(y_{i})}{u^{\prime }(y_{i})^{3}}(\ln (\unicode[STIX]{x1D709}_{i+1}-y_{i})-\ln (\unicode[STIX]{x1D709}_{i}-y_{i}))+o(1)\end{eqnarray}$$
(A 11b ) $$\begin{eqnarray}\displaystyle & \underset{c\rightarrow 0}{=} & \displaystyle \text{i}\unicode[STIX]{x03C0}\frac{u^{\prime \prime }(y_{i})}{u^{\prime }(y_{i})^{3}}+K_{R}+o(1),\end{eqnarray}$$
whereas we obtain
(A 12) $$\begin{eqnarray}\displaystyle K_{0,c}\underset{c\rightarrow 0}{=}-\frac{u^{\prime \prime }(0)}{u^{\prime }(0)^{3}}\left(-\!\ln \left(-\frac{c}{u^{\prime }(0)}\right)\right)+K_{R}+o(1). & & \displaystyle\end{eqnarray}$$

Finally, combining (A 6b ), (A 9), (A 11b ) and (A 12), we get the desired estimate of  $I_{c}$ :

(A 13) $$\begin{eqnarray}\displaystyle I_{c}\underset{c\rightarrow 0}{=}\frac{1}{u^{\prime }(0)}-c\frac{u^{\prime \prime }(0)}{u^{\prime }(0)^{3}}\ln \left(-\frac{c}{u^{\prime }(0)}\right)-c\text{i}\unicode[STIX]{x03C0}\mathop{\sum }_{i=1}^{m-1}\frac{u^{\prime \prime }(y_{i})}{u^{\prime }(y_{i})^{3}}+cK_{R}+o(c). & & \displaystyle\end{eqnarray}$$

If $u^{\prime }(0)=0$ , we cannot apply Vitali’s theorem to (A 1) for $i=0$ because the second and third terms diverge when $c\rightarrow 0$ . Hughes & Reid (Reference Hughes and Reid1965) computed the asymptotic expansion of $I_{c}$ for a special form of $u$ with $u^{\prime }(0)=0$ , but we rederive it using a different method for completeness. Letting $v=\sqrt{u}$ , we may write

(A 14) $$\begin{eqnarray}\displaystyle (u-c)^{2}=(v-\sqrt{c})(v+\sqrt{c}), & & \displaystyle\end{eqnarray}$$

with

(A 15a,b ) $$\begin{eqnarray}\displaystyle v^{\prime }(z)^{2}\xrightarrow[{z\rightarrow 0}]{}\frac{u^{\prime \prime }(0)}{2},\quad v^{\prime \prime }(z)\xrightarrow[{z\rightarrow 0}]{}\pm \sqrt{\frac{2}{u^{\prime \prime }(0)}}u^{\prime \prime \prime }(0), & & \displaystyle\end{eqnarray}$$

so that the above results can be applied to $v$ with $z^{+}$ and $z^{-}$ defined by $v(z^{+}(c))=\sqrt{c}$ and $v(z^{-}(c))=\sqrt{c}$ , leading to

(A 16) $$\begin{eqnarray}\displaystyle g^{\pm }(y,c)=\frac{1}{(v(y)\pm \sqrt{c})^{2}}-\frac{1}{v^{\prime }(z^{\pm }(c))^{2}(y-z^{\pm }(c))^{2}}+\frac{v^{\prime \prime }(z^{\pm }(c))}{v^{\prime }(z^{\pm }(c))^{3}(y-z^{\pm }(c))}. & & \displaystyle\end{eqnarray}$$

Since the functions $g^{\pm }$ are bounded by a constant independent of $c$ , we can now safely apply the Vitali theorem to the product $g^{+}g^{-}$ . Owing to the order of the different terms, it is sufficient to keep only one of them,

(A 17) $$\begin{eqnarray}\displaystyle I_{c}\underset{c\rightarrow 0}{{\sim}}-c\int _{0}^{\infty }\frac{1}{v^{\prime }(z^{+}(c))^{2}v^{\prime }(z^{-}(c))^{2}(y-z^{+}(c))^{2}(y-z^{-}(c))^{2}}, & & \displaystyle\end{eqnarray}$$

and after computing the residuals we finally obtain

(A 18) $$\begin{eqnarray}\displaystyle I_{c}\underset{c\rightarrow 0}{=}(8u^{\prime \prime }(0)c)^{-1/2}\text{i}\unicode[STIX]{x03C0}+o(c^{-1/2}). & & \displaystyle\end{eqnarray}$$

Appendix B. Validation of the solvers

Although the discretization methods used for this paper are relatively classical, the way the boundary conditions are imposed is new and it was thus necessary to conduct validation runs, which are reported here.

B.1 Navier–Stokes solver

As a test case for the Navier–Stokes solver, the set-up designed by Kramer et al. (Reference Kramer, Clercx and van Heijst2007) was considered. Contrary to the quadrupole set-up which has been studied in the body of the present paper, the dipole is not symmetric with respect to the channel centreline. The full span of the channel and the walls on both sides therefore needs to be taken into account. Two runs with $Re=650$ and $Re=2500$ were performed, respectively, with $512\times 255$ and $1024\times 511$ uniformly distributed grid points.

Table 3. Location $(x_{d},y_{d})$ and value $(\unicode[STIX]{x1D714}_{d})$ of maximum vorticity for various times and Reynolds numbers, compared to the data in table III of Kramer et al. (Reference Kramer, Clercx and van Heijst2007), for the dipole–wall test case.

In order to make a quantitative comparison, the same procedure as used by Kramer et al. (Reference Kramer, Clercx and van Heijst2007) is repeated here, namely to compare the location and amplitude of the main vortex core at several instants. The data are presented in table 3. Note that, at $t=0$ , there is a minor mistake in the reference data, since $x_{d}=0.1$ corresponds by construction to the location of the vorticity maximum for an isolated monopole, whereas the maximum of the dipole is slightly shifted due to interaction with the opposite-sign vortex. The results are otherwise in good agreement. The slight mismatch of the vortex position for $Re=2500$ might be due to the under-resolution of the present computation; nevertheless, the amplitude of the main vortex core agrees well. Detailed benchmarking is an interesting perspective for future studies, as many data for comparison are available in the literature (Clercx & Bruneau Reference Clercx and Bruneau2006; Kramer et al. Reference Kramer, Clercx and van Heijst2007).

B.2 Prandtl–Euler solver

For the Prandtl solver, the classical impulsively started cylinder studied by Van Dommelen & Shen (Reference Van Dommelen and Shen1980), henceforth VDS, in the Lagrangian framework is employed as a test case. It corresponds to a constant boundary pressure gradient given by

(B 1) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x2202}_{x}p_{E}(x,0,t)=-\!\sin (x)\cos (x), & & \displaystyle\end{eqnarray}$$

and an initial vorticity which is a Dirac distribution,

(B 2) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D714}_{P}(x,y,0)=-\unicode[STIX]{x1D6FF}_{0}(y)\sin (x). & & \displaystyle\end{eqnarray}$$

For spatial discretization, $1023$ grid points are considered on the interval $[0,\unicode[STIX]{x03C0}]$ in the $x$ direction, taking advantage of the odd symmetry of the solution, and $513$ grid points with $L_{y}=32$ in the $y$ direction. The initial Dirac distribution is approximated by letting

(B 3a,b ) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D714}_{P}(x_{i},0,0)=-\frac{\sin (x_{i})}{\unicode[STIX]{x0394}x},\quad \unicode[STIX]{x1D714}_{P}(x_{i},y,0)=0\quad \text{if }y>0. & & \displaystyle\end{eqnarray}$$

The results are then compared with figure 10 and table II of VDS. Figure 15 shows the wall shear stress for different time instants. It is in very good qualitative agreement with figure 10 of VDS, except maybe at very short times, which is to be expected given the singular initial condition. Additionally, our estimate for the quantity $F_{w}^{\prime \prime }$ at $t=3$ is $1.1061$ , which is in good agreement with the value $1.1122$ found by VDS at their much lower resolution (in doing this comparison we have assumed that the undefined quantity $T$ appearing in table II of VDS corresponds to $t/2$ ).

Figure 15. Wall shear stress at different instants for the impulsively started cylinder test case.

B.3 Orr–Sommerfeld solver

To validate the MATLAB code used to compute Orr–Sommerfeld eigenvalues, we use the Blasius boundary layer as a test case. The obtained marginal stability curve (figure 16) is in good agreement with figure 16.11 provided by Schlichting (Reference Schlichting1979).

Figure 16. Unstable region for the Blasius boundary layer.

Figure 17. Vorticity fields for the Navier–Stokes solution (left) and the Euler solution (right) in the bulk flow at $t=30,50$ and $55.1$ . The corresponding vorticity fields in the boundary layer for the Navier–Stokes (left) and the Prandtl solutions (right) are given directly below. The corresponding movie can be downloaded from https://doi.org/10.1017/jfm.2018.396.

Figure 18. Continuation of figure 17 at time instants $t=57$ and $100$ . The Prandtl solution is singular and thus not shown. The corresponding movie can be downloaded from https://doi.org/10.1017/jfm.2018.396.

Appendix C. Overall flow evolution: comparison of Navier–Stokes and Euler/Prandtl flows

The dipole first shoots towards the lower channel wall. The Navier–Stokes vorticity field in the bulk (figure 17 top, left) remains very close to the Euler vorticity field (figure 17 top, right). Plotting the Navier–Stokes flow in the Prandtl boundary layer units (figure 17 top, left) reveals that it is smooth, and very well approximated by the solution of the Prandtl equations (figure 17 top, right). The vorticity along the boundary (figure 19 a) converges to the Prandtl values as the Reynolds number increases.

As the dipole approaches the wall (figure 17, middle), the pressure gradient along the wall becomes more intense and steeper, which causes a strong inward diffusion of vorticity at the wall, as well as increased vorticity gradients within the boundary layer. At time $t=50$ , we still observe convergence of the Navier–Stokes solution at high Reynolds number towards the Euler flow in the bulk and towards the Prandtl flow in the boundary layer. However, looking at the vorticity along the boundary (figure 19 b) reveals a larger difference between Prandtl and Navier–Stokes flows than at $t=30$ .

As the Prandtl solution approaches its singularity time $t^{\ast }\approx 55.6$ (figure 17, bottom), parallel vorticity gradients increase rapidly, and soon the cut-off parallel wavenumber of the numerical scheme becomes insufficient to resolve it. The convergence of the Navier–Stokes boundary vorticity is lost over a wide interval in $x$ , and the vorticity around $x=0.61$ adopts a stronger scaling with Reynolds number (figure 19 c).

After the singularity, at $t=57$ (figure 18, top), oscillations in the vorticity appear, while the bulk flow still looks similar for Navier–Stokes and Euler. Following the new vorticity extremum that has appeared at the boundary, a cascade of extrema with opposite signs appear (for sufficiently high Reynolds number), exciting increasingly fine scales parallel to the wall (figure 19 d).

Figure 19. Vorticity along the boundary at $t=30,~50,~55.1$ and $57$ .

At much later times, $t=100$ (figure 18, bottom), the Euler and Navier–Stokes solutions have become completely different. In the Euler case, the vortices glide along the wall, having paired with their mirror image, no new vorticity has been generated and the energy is conserved. In the Navier–Stokes case, the detachment process has led to the formation of two new vortices (shown in cyan and magenta in figure 18, bottom, left) of much larger amplitude than those of the incoming dipole. The activity in the boundary layer remains intense, leading to the ejection of smaller structures.

References

Balay, S., Adams, M. F., Brown, J., Brune, P., Buschelman, K., Eijkhout, V., Gropp, W. D., Kaushik, D., Knepley, M. G., McInnes, L. C., Rupp, K., Smith, B. F. & Zhang, H.2013 PETSc users manual. Tech. Rep. ANL-95/11, Revision 3.4. Argonne National Laboratory.Google Scholar
Bardos, C. & Titi, E. S. 2013 Mathematics and turbulence: where do we stand? J. Turbul. 14 (3), 4276.Google Scholar
Bedrossian, J., Masmoudi, N. & Vicol, V. 2016 Enhanced dissipation and inviscid damping in the inviscid limit of the Navier–Stokes equations near the two dimensional Couette flow. Arch. Rat. Mech. Anal. 219 (3), 10871159.Google Scholar
Blasius, H. 1908 Grenzschichten in Flüssigkeiten mit kleiner Reibung. Z. Math. Phys. 56, 137; Engl. transl. in NACA TM 1256 (1950).Google Scholar
Bowles, R. 2006 Lighthill and the triple-deck, separation and transition. J. Engng Maths 56 (4), 445460.Google Scholar
Bowles, R. I., Davies, C. & Smith, F. T. 2003 On the spiking stages in deep transition and unsteady separation. J. Engng Maths 45 (3), 227245.Google Scholar
Brinckman, K. W. & Walker, J. D. A. 2001 Instability in a viscous flow driven by streamwise vortices. J. Fluid Mech. 432, 127166.Google Scholar
Cassel, K. W. 2000 A comparison of Navier–Stokes solutions with the theoretical description of unsteady separation. Phil. Trans. R. Soc. Lond. A 358 (1777), 32073227.Google Scholar
Cassel, K. W. & Obabko, A. V. 2010 A Rayleigh instability in a vortex-induced unsteady boundary layer. Phys. Scr. T142, 114.Google Scholar
Clercx, H. & Bruneau, G.-H. 2006 The normal and oblique collision of a dipole with a no-slip boundary. Comput. Fluids 35, 245279.Google Scholar
Clercx, H. & van Heijst, G. J. 2002 Dissipation of kinetic energy in two-dimensional bounded flows. Phys. Rev. E 65, 066305.Google Scholar
Clercx, H. & van Heijst, G. J. 2017 Dissipation of coherent structures in confined two-dimensional turbulence. Phys. Fluids 29, 111103.Google Scholar
Cowley, S. J. 2002 Laminar boundary-layer theory: a 20th century paradox? In Mechanics for a New Millennium (ed. Aref, H. & Phillips, J. W.), pp. 389412. Kluwer.Google Scholar
Demmel, J. W., Gilbert, J. R. & Li, X. S. 1999 An asynchronous parallel supernodal algorithm for sparse Gaussian elimination. SIAM J. Matrix Anal. Applics 20 (4), 915952.Google Scholar
E, W. & Engquist, B. 1997 Blowup of solutions of the unsteady Prandtl’s equation. Commun. Pure Appl. Maths 50 (12), 12871293.Google Scholar
E, W. & Liu, J. G. 1996 Vorticity boundary condition and related issues for finite difference schemes. J. Comput. Phys. 124 (2), 368382.Google Scholar
Elliott, J. W., Smith, F. T. & Cowley, S. J. 1983 Breakdown of boundary layers: (i) on moving surfaces; (ii) in semi-similar unsteady flow; (iii) in fully unsteady flow. Geophys. Astrophys. Fluid Dyn. 25 (1–2), 77138.Google Scholar
Gamet, L., Ducros, F., Nicoud, F. & Poinsot, T. 1999 Compact finite difference schemes on non-uniform meshes: application to direct numerical simulations of compressible flows. Intl J. Numer. Meth. Fluids 29 (2), 159191.Google Scholar
Gargano, F., Sammartino, M. & Sciacca, V. 2009 Singularity formation for Prandtl’s equations. Physica D 238 (19), 19751991.Google Scholar
Gargano, F., Sammartino, M. & Sciacca, V. 2011 High Reynolds number Navier–Stokes solutions and boundary layer separation induced by a rectilinear vortex. Comput. Fluids 52, 7391.Google Scholar
Gargano, F., Sammartino, M., Sciacca, V. & Cassel, K. W. 2014 Analysis of complex singularities in high-Reynolds-number Navier–Stokes solutions. J. Fluid Mech. 747, 381421.Google Scholar
Gérard-Varet, D. & Dormy, E. 2010 On the ill-posedness of the Prandtl equation. J. Am. Math. Soc. 23 (2), 591609.Google Scholar
Gérard-Varet, D. & Masmoudi, N. 2015 Well-posedness for the Prandtl system without analyticity or monotonicity. Ann. Sci. Éc. Norm. Supér. 48 (6), 12731325.Google Scholar
Goldstein, S. 1948 On laminar boundary-layer flow near a position of separation. Q. J. Mech. Appl. Math. 1 (1), 4369.Google Scholar
Grenier, E. 2000 On the nonlinear instability of Euler and Prandtl equations. Commun. Pure Appl. Math. 53, 10671091.Google Scholar
Grenier, E., Guo, Y. & Nguyen, T. T. 2015 Spectral stability of Prandtl boundary layers: an overview. Analysis 35 (4), 343355.Google Scholar
Grenier, E., Guo, Y. & Nguyen, T. T. 2016 Spectral instability of characteristic boundary layer flows. Duke Math. J. 165 (16), 30853146.Google Scholar
Grenier, E. & Nguyen, T. T.2018 Sublayer of Prandtl boundary layers. Arch. Rat. Mech. Anal. doi:10.1007/s00205-018-1235-3.Google Scholar
Gresho, P. M. 1991 Incompressible fluid dynamics: some fundamental formulation issues. Annu. Rev. Fluid Mech. 23 (1), 413453.Google Scholar
Heisenberg, W. 1924 Über Stabilität und Turbulenz von Flüssigkeitsströmen. Ann. Physik 74, 577627.Google Scholar
Hughes, T. H. & Reid, W. H. 1965 The stability of laminar boundary layers at separation. J. Fluid Mech. 23 (04), 737747.Google Scholar
Ingham, D. B. 1984 Unsteady separation. J. Comput. Phys. 53 (1), 9099.Google Scholar
von Kármán, T. 1921 Über laminare und turbulente Reibung. Z. Angew. Math. Mech. 1 (4), 233252; Engl. transl. in NACA TM 1092 (1946).Google Scholar
Kato, T. 1972 Nonstationary flows of viscous and ideal fluids in ℝ3 . J. Funct. Anal. 9, 296305.Google Scholar
Kato, T. 1984 Remarks on zero viscosity limit for nonstationary Navier–Stokes flows with boundary. In Seminar on Nonlinear Partial Differential Equations, pp. 8598. MSRI.Google Scholar
Keetels, G. H., Kramer, W., Clercx, H. & van Heijst, G. J. 2011 On the Reynolds number scaling of vorticity production at no-slip walls during vortex–wall collisions. Theor. Comput. Fluid Dyn. 25, 293300.Google Scholar
Kelliher, J. P. 2007 On Kato’s conditions for vanishing viscosity. Indiana Univ. Math. J. 56 (4), 17111722.Google Scholar
Kramer, W., Clercx, H. & van Heijst, G. J. 2007 Vorticity dynamics of a dipole colliding with a no-slip wall. Phys. Fluids 19, 126603.Google Scholar
Lele, S. K. 1992 Compact finite difference schemes with spectral-like resolution. J. Comput. Phys. 103 (1), 1642.Google Scholar
Li, X. S., Demmel, J. W., Gilbert, J. R., Grigori, L., Shao, M. & Yamazaki, I.1999 SuperLU Users’ Guide. Tech. Rep. LBNL-44289. Lawrence Berkeley National Laboratory, http://crd-legacy.lbl.gov/∼xiaoye/SuperLU/superlu_ug.pdf. Last update: August 2011.Google Scholar
Lin, C. C. 1967 The Theory of Hydrodynamic Stability. Cambridge University Press.Google Scholar
Maekawa, Y. 2013 Solution formula for the vorticity equations in the half-plane with application to high vorticity creation at zero viscosity limit. Adv. Differ. Equ. 18 (1/2), 101146.Google Scholar
Maekawa, Y. 2014 On the inviscid limit problem of the vorticity equations for viscous incompressible flows in the half-plane. Commun. Pure Appl. Maths 67 (7), 10451128.Google Scholar
Marusic, I., McKeon, B. J., Monkewitz, P. A., Nagib, H. M., Smits, A. J. & Sreenivasan, K. R. 2010 Wall-bounded turbulent flows at high Reynolds numbers: recent advances and key issues. Phys. Fluids 22 (6), 065103.Google Scholar
Nguyen van yen, R., Farge, M. & Schneider, K. 2009 Wavelet regularization of a Fourier–Galerkin method for solving the 2D incompressible Euler equations. ESAIM: Proc. 29, 89107.Google Scholar
Nguyen van yen, R., Farge, M. & Schneider, K. 2011 Energy dissipating structures produced by walls in two-dimensional flows at vanishing viscosity. Phys. Rev. Lett. 106 (18), 184502.Google Scholar
Nguyen van yen, R., Kolomenskiy, D. & Schneider, K. 2014 Approximation of the Laplace and Stokes operators with Dirichlet boundary conditions through volume penalization: a spectral viewpoint. Numer. Math. 128, 301338.Google Scholar
Obabko, A. V. & Cassel, K. W. 2002 Navier–Stokes solutions of unsteady separation induced by a vortex. J. Fluid Mech. 465, 99130.Google Scholar
Oleinik, O. A. 1966 On the mathematical theory of boundary layer for an unsteady flow of incompressible fluid. J. Appl. Math. Mech. 30 (5), 951974.Google Scholar
Orlandi, P. 1990 Vortex dipole rebound from a wall. Phys. Fluids A 2, 14291436.Google Scholar
Peridier, V. J., Smith, F. T. & Walker, J. D. A. 1991a Vortex-induced boundary-layer separation. Part 1. The unsteady limit problem Re . J. Fluid Mech. 232, 99131.Google Scholar
Peridier, V. J., Smith, F. T. & Walker, J. D. A. 1991b Vortex-induced boundary-layer separation. Part 2. Unsteady interacting boundary-layer theory. J. Fluid Mech. 232, 133165.Google Scholar
Prandtl, L. 1904 Über Flüssigkeitsbewegung bei sehr kleiner Reibung. In Proceedings of the 3rd International Mathematical Congress, Heidelberg, pp. 484491.Google Scholar
Prandtl, L. 1921 Bemerkungen über die Entstehung der Turbulenz. Z. Angew. Math. Mech. 1 (6), 431436.Google Scholar
Rayleigh, Lord 1880 On the stability, or instability, of certain fluid motions. Proc. Lond. Math. Soc. 11 (1), 5770.Google Scholar
Reed, H. L., Saric, W. S. & Arnal, D. 1996 Linear stability theory applied to boundary layers. Annu. Rev. Fluid Mech. 28 (1), 389428.Google Scholar
Robinson, S. K. 1991 Coherent motions in the turbulent boundary layer. Annu. Rev. Fluid Mech. 23, 601639.Google Scholar
le Rond d’Alembert, J.-B. 1768 Paradoxe proposé aux géomètres sur la résistance des fluides. In Opuscules Mathématiques, vol. 5, chap. XXXIV, pp. 132138. Briasson.Google Scholar
Sammartino, M. & Caflisch, R. E. 1998 Zero viscosity limit for analytic solutions of the Navier–Stokes equation on a half-space. Commun. Math. Phys. 192, 433491.Google Scholar
Schlichting, H. 1979 Boundary Layer Theory. McGraw-Hill.Google Scholar
Schubauer, G. B. & Skramstad, H. K. 1947 Laminar boundary layer oscillations and transitions on a flat plate. J. Aero. Sci. 14, 6976; also NACA Report 909 (1948).Google Scholar
Sears, W. R. & Telionis, D. P. 1975 Boundary-layer separation in unsteady flow. SIAM J. Appl. Maths 28 (1), 215235.Google Scholar
Smith, F. T. 1988 Finite-time break-up can occur in any unsteady interacting boundary layer. Mathematika 35 (02), 256273.Google Scholar
Sommerfeld, A. 1909 Ein Beitrag zur hydrodynamischen Erklärung der turbulenten Flüssigkeitsbewegung. In Proceedings of the 4th International Mathematical Congress, Rome 1908, vol. 3, pp. 116124.Google Scholar
Stewartson, K. 1974 Multistructured boundary layers on flat plates and related bodies. In Adv. Appl. Mech., vol. 14, pp. 145239. Academic Press.Google Scholar
Sutherland, D., Macaskill, C. & Dritschel, D. G. 2013 The effect of slip length on vortex rebound from a rigid boundary. Phys. Fluids 25 (9), 093104.Google Scholar
Temam, R. & Wang, X. 1997 The convergence of the solutions of the Navier–Stokes equations to that of the Euler equations. Appl. Math. Lett. 10 (5), 2933.Google Scholar
Tietjens, O. 1925 Beiträge zur Entstehung der Turbulenz. Z. Angew. Math. Mech. 5 (3), 200217.Google Scholar
Tollmien, W. 1929 Über die Entstehung der Turbulenz. Nachr. Ges. Wiss. Göttingen, Math.-Physik. Klasse 2144; Engl. transl. in NACA Tech. Rep. 609 (1931).Google Scholar
Van Dommelen, L.1981 Unsteady boundary layer separation. PhD thesis, Cornell University.Google Scholar
Van Dommelen, L. L. & Cowley, S. J. 1990 On the Lagrangian description of unsteady boundary-layer separation. Part 1. General theory. J. Fluid Mech. 210, 593626.Google Scholar
Van Dommelen, L. L. & Shen, S. F. 1980 The spontaneous generation of the singularity in a separating laminar boundary layer. J. Comput. Phys. 38 (2), 125140.Google Scholar
Figure 0

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.

Figure 1

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

Figure 2

Table 1. Definitions and initial values of several integral quantities of interest.

Figure 3

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

Figure 4

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.

Figure 5

Figure 4. Comparison between Navier–Stokes at $Re=123\,075$ (black) and Euler/Prandtl (blue). (a,b) Contour lines of vorticity field at $t=30$ (a) and $t=50$ (b), for $\unicode[STIX]{x1D714}=0.1,0.2,\ldots ,1$. (c,d) Contour lines of corresponding boundary layer vorticity fields in Prandtl variables, for $\unicode[STIX]{x1D714}=-|\unicode[STIX]{x1D714}|_{max}/6,-2|\unicode[STIX]{x1D714}|_{max}/6,\ldots ,-5|\unicode[STIX]{x1D714}|_{max}/6$, where $|\unicode[STIX]{x1D714}|_{max}=3.75\times 10^{-4}$ and $|\unicode[STIX]{x1D714}|_{max}=8.70\times 10^{-3}$, respectively.

Figure 6

Figure 5. Contour lines of vorticity field in Prandtl variables at $t=54$ (a) and $t=55.3$ (b), for $\unicode[STIX]{x1D714}=5|\unicode[STIX]{x1D714}|_{max}/6,\ldots ,|\unicode[STIX]{x1D714}|_{max}/6,0,-|\unicode[STIX]{x1D714}|_{max}/6,-2|\unicode[STIX]{x1D714}|_{max}/6,\ldots ,-5|\unicode[STIX]{x1D714}|_{max}/6$, where $|\unicode[STIX]{x1D714}|_{max}=1.25\times 10^{-2}$ and $|\unicode[STIX]{x1D714}|_{max}=1.30\times 10^{-2}$, respectively. Black contour lines correspond to Navier–Stokes at $Re=123\,075$, while blue lines correspond to Euler/Prandtl. The $\unicode[STIX]{x1D714}=0$ contour line is shown in bold.

Figure 7

Figure 6. Comparison between Navier–Stokes (black), Prandtl–Euler (blue) and interactive boundary layer (red) models at $t=50$ (a) and $t=54$ (b). Contour lines of vorticity (a,b) and vorticity along the boundary (c,d). The $\unicode[STIX]{x1D714}=0$ contour line is shown in bold.

Figure 8

Figure 7. Vorticity along the boundary for Navier–Stokes (solid line) for varying Reynolds number (from $Re=7692$ to $123\,075$) and Prandtl (dashed line).

Figure 9

Figure 8. Fourier spectra of boundary vorticity for Navier–Stokes (solid line) at different Reynolds numbers (from $Re=7692$ to $123\,075$) and Prandtl (dashed line). The spectra have been smoothed to eliminate fast oscillations due to the odd symmetry of the function (see text).

Figure 10

Figure 9. Stability boundary subject to the frozen flow approximation, shown as a function of $x$ at $t=54$ and for $Re=123\,075$. The upper and lower marginal wavenumbers estimated from our asymptotic analysis (2.51) and (2.54) are shown in dashed lines.

Figure 11

Figure 10. Contour lines of vorticity field in Prandtl variables for Navier–Stokes at $t=56$ (a) and $t=56.9$ (b). The $\unicode[STIX]{x1D714}=0$ contour line is shown in bold.

Figure 12

Figure 11. Vorticity along the boundary for various Reynolds numbers (from $Re=7692$ to $123\,075$) in Prandtl variables at $t=56$ (a) and $t=56.9$ (b).

Figure 13

Figure 12. Reynolds-number scalings of maximum vorticity (a) and enstrophy (b) for Prandtl and Navier–Stokes before detachment ($t=53$), and for Navier–Stokes after detachment ($t=56.9$).

Figure 14

Figure 13. Contour lines of vorticity field at $t=80$ for Navier–Stokes at $Re=15\,384$ and Euler/Prandtl solutions.

Figure 15

Figure 14. Comparison between Navier–Stokes solutions at the same Reynolds number ($Re=123\,075$) with two numerical resolutions, in order to check for convergence: (a) $t=56$, (b) $t=56.9$.

Figure 16

Table 3. Location $(x_{d},y_{d})$ and value $(\unicode[STIX]{x1D714}_{d})$ of maximum vorticity for various times and Reynolds numbers, compared to the data in table III of Kramer et al. (2007), for the dipole–wall test case.

Figure 17

Figure 15. Wall shear stress at different instants for the impulsively started cylinder test case.

Figure 18

Figure 16. Unstable region for the Blasius boundary layer.

Figure 19

Figure 17. Vorticity fields for the Navier–Stokes solution (left) and the Euler solution (right) in the bulk flow at $t=30,50$ and $55.1$. The corresponding vorticity fields in the boundary layer for the Navier–Stokes (left) and the Prandtl solutions (right) are given directly below. The corresponding movie can be downloaded from https://doi.org/10.1017/jfm.2018.396.

Figure 20

Figure 18. Continuation of figure 17 at time instants $t=57$ and $100$. The Prandtl solution is singular and thus not shown. The corresponding movie can be downloaded from https://doi.org/10.1017/jfm.2018.396.

Figure 21

Figure 19. Vorticity along the boundary at $t=30,~50,~55.1$ and $57$.

Nguyen van yen et al. supplementary movie

Time evolution, from $t=0$ to $90$, of the vorticity field in the bulk flow for the Euler solution (left) and the Navier--Stokes solution (right). The Prandtl solution becomes singular at $t^{\star}=55.6$ and is thereafter not shown. The corresponding vorticity field in the boundary layer for the Prandtl solution (left) and the Navier--Stokes solution (right) are given directly below.

Download Nguyen van yen et al. supplementary movie(Video)
Video 5 MB