#  Energy dissipation caused
by boundary layer instability
at vanishing viscosity

Romain Nguyen van yen1 Mathias Waidmann1 Rupert Klein1 Marie Farge2 Kai
Schneider3 1 Institut für Mathematik, Freie Universität Berlin, Arnimallee 6,
14185 Berlin, Germany
2 LMD-CNRS, Ecole Normale Supérieure, 24 rue Lhomond, 75231 Paris Cedex 5,
France
3 Institut de Mathématiques de Marseille, Aix-Marseille Université & CNRS,
Marseille, France

###### 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 $\mathrm{Re}^{3/8}$ and $\mathrm{Re}^{1/2}$,
respectively. By linear superposition the unstable modes induce a boundary
vorticity flux of order $\mathrm{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.

###### Contents

  1. 1 Introduction
  2. 2 Models
    1. 2.1 Navier–Stokes model
    2. 2.2 Prandtl–Euler model
    3. 2.3 Interactive boundary layer model
    4. 2.4 Orr–Sommerfeld model
      1. 2.4.1 Formulation
      2. 2.4.2 Inviscid mode
      3. 2.4.3 Viscous correction
      4. 2.4.4 Construction of unstable modes
      5. 2.4.5 The case $u_{P}^{\prime}(0)=0$
      6. 2.4.6 Physical interpretation
  3. 3 Solvers
    1. 3.1 Setup
    2. 3.2 Navier–Stokes solver
    3. 3.3 Prandtl–Euler solver
    4. 3.4 Interactive boundary layer solver
    5. 3.5 Orr–Sommerfeld solver
  4. 4 Results
    1. 4.1 Before detachment
    2. 4.2 Prandtl blow-up
    3. 4.3 Large scale interaction and instability
    4. 4.4 Detachment and production of dissipative structures
    5. 4.5 Later evolution
    6. 4.6 Convergence checks
  5. 5 Discussion
  6. A Estimation of the integral $\boldmath I_{c}$
  7. B Validation of the solvers
    1. B.1 Navier–Stokes solver
    2. B.2 Prandtl–Euler solver
    3. B.3 Orr–Sommerfeld solver
  8. C Overall flow evolution: comparison of Navier–Stokes and Euler/Prandtl flows

##  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 is one of the central unknowns in the field of fluid mechanics.
From the vorticity transport equations he had derived, d’Alembert (1768)
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 19th
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 $\nu$ of the fluid,
or equivalently, after non-dimensionalization, on the inverse of the Reynolds
number $\mathrm{Re}$, predict that friction should become negligible when
$\mathrm{Re}\gg 1$.

The group working at Göttingen University, notably Prandtl (1904) and Blasius
(1908), were first to come up with a way of computing, under some hypotheses,
the asymptotic behavior of the force in the limit $\mathrm{Re}\to\infty$.
Their method relies on the notion that viscous effects are confined to a layer
of thickness of order $\mathrm{Re}^{-1/2}$ along the boundary, now called
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(\mathrm{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 $\mathrm{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, 1979).

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 neighborhood of the wall
reverses direction at this point. After the formal asymptotic expansion
achieved by Goldstein (1948), which clearly indicates a singularity, a viscous
scaling of the parallel coordinate was developed to analyze the flow near the
separation point, finally leading to the so-called triple deck structure
(Stewartson, 1974). This steady asymptotic theory is instructive and well
developed but remains unsatisfactory, since in practice all flows are unsteady
above a certain Reynolds number.

The classical understanding of the onset of unsteadiness in shear flows goes
back to Rayleigh (1880), who showed in particular that a necessary condition
for inviscid instability is that the base velocity profile has an inflexion
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 (1909) 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
nonzero second derivative of the base flow could trigger an instability which
is absent in the inviscid setting (Prandtl, 1921; Tietjens, 1925). Following
some of the ideas developed by Heisenberg (1924) in his thesis, Tollmien
(1929) 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.

Beyond this, one is faced with the full unsteady Navier–Stokes equations (NSE)
at high Reynolds number, and the BL problem becomes so wide reaching that it
has been investigated almost independently for the past 60 years by two
disctinct 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 lead in the 1950s to the establishment of
the Moore-Rott-Sears criterion (Sears & Telionis, 1975), 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
(1975) conjectured that separation coincides with the appearance of a
singularity in the solution at some finite time $t^{*}$. The approach to such
a singularity was confirmed using numerical experiments by Van Dommelen & Shen
(1980) for the impulsively started cylinder, and then analyzed in detail by
power series expansions in Lagrangian variables by Van Dommelen & Cowley
(1990). We refer the reader to van Dommelen’s PhD thesis (1981) for a detailed
pedagogical account of these findings. Later work largely supports his results
(Ingham, 1984; E & Engquist, 1997; Gargano et al., 2009).

The next challenge was to understand what happens to the actual NSE flow as
the corresponding Prandtl solution becomes singular. Elliott et al. (1983)
obtained the estimate $(t^{*}-t)=O(\mathrm{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, 1988; Peridier et al., 1991a, 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 on the vanishing viscosity limit for solutions
of the NSE on the other hand. Local well-posedness for the Prandtl equations
was proved long ago by Oleinik (1966) in cases where detachment is not
expected, i.e., for monotonous initial data and favorable pressure gradient.
Sammartino & Caflisch (1998) 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 Gerard-Varet & Masmoudi (2013), but the required decay of
Fourier coefficients is still faster than any power of $k$. Maekawa (2014)
proved the convergence of the Euler and Prandtl solution versus the
Navier–Stokes solution in the $L^{\infty}$ norm with order $\sqrt{\nu}$ for an
initial condition where the vorticity is compactly supported at finite
distance from the wall.

Kato (1984) made a decisive contribution, by linking the vanishing viscosity
limit problem to the behavior of energy dissipation at the boundary. Kato’s
theorem implies, in the particular case of a flow in a smooth two-dimensional
domain $\Omega$ with smooth initial data and without forcing, the equivalence
between the following assertions:

  1. 1.

the NS flow converges to the Euler flow uniformly in time in the energy norm,

  2. 2.

the energy dissipation associated to the NS flow, integrated over a strip of
thickness proportional to $\mathrm{Re}^{-1}$ around the solid, which we will
call the Kato layer, tends to zero as $\mathrm{Re}\to\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 $\mathrm{Re}^{-1}$ for detachment to
be possible. Later refinements of Kato’s work linked breakdown to scalings
$O(\mathrm{Re}^{\alpha})$ with $\alpha\geq 3/4$ of the wall pressure gradient
(Temam & Wang, 1997), and $O(\mathrm{Re}^{\beta}$ with $\beta\geq 1/2$ of the
$L^{2}$ norm of velocity in the Kato layer (Kelliher, 2007).

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 (1990), 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 underly 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 $\mathrm{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
(2000) 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 (2010) 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. Once again, the details of the proof show that the ill-posedness
is due to modes with large wavenumber in the direction parallel to the wall.
More recently, Grenier et al. (2014b, a, 2015) started working on the
ambitious program aiming to achieve a rigorous mathematical description of
instabilities in generic BL flows.

Several fluid dynamicists also advanced the idea that instability-type
mechanisms may play an important role in the process of unsteady detachment.
Cassel (2000) directly compared numerical solutions of the NSE 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 (1980), developing strong parallel
velocity gradients in a process which seems to lead to a finite time
singularity associated to 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 behavior, characterized by
the appearance of strong oscillations in the wall-parallel pressure gradient,
which he was not able to explain. Brinckman & Walker (2001) also saw
oscillations, and claimed that they were due to a Rayleigh instability of the
shear velocity profile, an hypothesis which was developed further in the
review paper by Cowley (2002). Although the numerics underlying these
findings, and therefore their interpretation in terms of a Rayleigh
instability, were later invalidated due to their insufficient grid resolution
by Obabko & Cassel (2002), the existence of an instability mechanism has
continued to be a hot topic in later papers on unsteady detachment (Bowles et
al., 2003; Bowles, 2006; Cassel & Obabko, 2010; Gargano et al., 2011, 2014).
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 to high Reynolds number flows, even in
2D, 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 et al., 2011), we computed a series of
dipole-wall collisions, a well studied academic flow introduced by Orlandi
(1990). 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 to provide 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 $\mathrm{Re}^{-1}$ (for a more detailed analysis
of the scheme, see Nguyen van yen et al., 2014). In this setting, we did find
indications that energy dissipation converges to a finite value when
$\mathrm{Re}\to\infty$.

Sutherland et al. (2013) then confirmed our findings using a Chebychev method
with exact Navier boundary conditions, but in the no-slip case concluded that
energy dissipation vanishes, in accordance with an earlier claim of Clercx &
van Heijst (2002). This is quite unexpetected since, due to the spectral
properties of the Stokes operator, the no-slip boundary condition is stiffer
that any Navier boundary condition with a nonzero slip length, and should thus
generate larger gradients or in other words more dissipation. In fact, looking
more closely at their results, it appears that the central claim is based on a
single computation (see Figs. 17 and 18 of their paper), for which a
convergence test is not provided, which in our opinion leaves the matter
unsettled.

This is why we propose 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 occurence of dissipation. We
then check that all these processes are actually ocurring in our numerical
solution of the dipole-wall initial value problem. For this, adopting a
methodology similar to the one employed by Cassel (2000), we proceed by direct
comparison of the NS solution and of the corresponding Prandtl approximation.

In the first section, 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 the second section, we introduce the model discretizations which we have
used for our numerical computations. In the third section we describe the
numerical results. Finally, we analyze 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 $\Omega$
read

$\displaystyle\partial_{t}\mathbf{u}+\left(\mathbf{u}\cdot\boldsymbol{\nabla}\right)\mathbf{u}=\nu\Delta\mathbf{u}-\boldsymbol{\nabla}p$
(1a) $\displaystyle\boldsymbol{\nabla}\cdot\mathbf{u}=0,$ (1b) where
$\mathbf{u}(\mathbf{x},t)$ is the velocity field, $p$ is the pressure field,
$\nu$ is the kinematic viscosity, and we shall denote $(u,v)$ the two
components of $\mathbf{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

$\Omega=\left\\{(x,y)\mid x\in\mathbb{T},y\in]0,1[\right\\},$ (1c)

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$,

$\mathbf{u}(x,0,t)=\mathbf{u}(x,1,t)=0,$ (1d)

and we specify an initial condition

$\mathbf{u}(x,y,0)=\mathbf{u}_{i}(x,y),$ (1e)

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:

$\mathrm{Re}=\frac{UL}{\nu}$ (2)

When discretizing this system, difficulties arise due to the interplay between
the divergence condition (1b) and the no-slip boundary condition (1d). We have
chosen to work with the vorticity formulation of the NSE, which eliminates the
divergence constraint at the cost of transforming the Dirichlet boundary
conditions on $u$ into a non-local integral constraint on $\omega$. Although
they used to be controversial, such formulations are now well established (see
Gresho, 1991; E & Liu, 1996; Maekawa, 2013), 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 $\omega=\partial_{x}v-\partial_{y}u$ satisfies the
transport equation

$\partial_{t}\omega+\boldsymbol{\nabla}\cdot(\mathbf{u}\omega)=\nu\Delta\omega,$
(3)

with initial data

$\omega_{i}(x,y):=\boldsymbol{\nabla}\times\mathbf{u}_{i}(x,y),$ (4)

where $\mathbf{u}$ is expressed as a function of $\omega$ by means of the
stream function $\psi$ defined by

$\displaystyle u$ $\displaystyle=-\partial_{y}\psi$ (5a) $\displaystyle v$
$\displaystyle=\partial_{x}\psi,$ (5b)

which in turn satisfies the Poisson equation

$\Delta\psi=\omega.$ (6)

From the wall-normal component of (1d) and the fact that
$\int_{\Omega}\mathbf{u}=0$ is a constant of motion, a Dirichlet boundary
condition for $\psi$ follows,

$\psi(x,0,t)=\psi(x,1,t)=0,$ (7)

which uniquely determines $\psi$, and therefore $\mathbf{u}$, as a function of
$\omega$.

To close the problem, the tangential component of (1d), which has not yet been
used, needs to be reformulated into the missing boundary condition on $\omega$
necessary because of the presence of a Laplacian in (3). A general discussion
of this issue has been carried out by Gresho (1991). In our case, due to the
simple geometry, (5)-(7) can be solved explicitly to get an expression for
$\psi$. For this, we first introduce the Fourier coefficients

$\widehat{\omega}_{k}(y)=\int_{0}^{1}\omega(x,y)e^{-\iota\kappa
x}\mathrm{d}x,$ (8)

where $\kappa=2\pi k$, and the corresponding reconstruction formula

$\omega(x,y)=\sum_{k\in\mathbb{Z}}\widehat{\omega}_{k}(y)e^{\iota\kappa x},$
(9)

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

$-\kappa^{2}\widehat{\psi}_{k}+\partial_{y}^{2}\widehat{\psi}_{k}=\widehat{\omega}_{k}.$
(10)

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

$\widehat{\psi}_{k}(y)=-\frac{1}{2|\kappa|\left(1-e^{-2|\kappa|}\right)}\Bigg{(}\left(1-e^{-2|\kappa|y}\right)\int_{y}^{1}\widehat{\omega}_{k}(y^{\prime})\left(e^{-|\kappa|(y^{\prime}-y)}-e^{-|\kappa|(2-y^{\prime}-y)}\right)\mathrm{d}y^{\prime}\\\
+\left(1-e^{-2|\kappa|(1-y)}\right)\int_{0}^{y}\widehat{\omega}_{k}(y^{\prime})\left(e^{-|\kappa|(y-y^{\prime})}-e^{-|\kappa|(y+y^{\prime})}\right)\mathrm{d}y^{\prime}\Bigg{)},$
(11a) for $k\neq 0$, and
$\widehat{\psi}_{0}(y)=y\int_{0}^{1}(1-y^{\prime})\widehat{\omega}_{0}(y^{\prime})\mathrm{d}y^{\prime}+\int_{0}^{y}(y-y^{\prime})\widehat{\omega}_{0}(y^{\prime})\mathrm{d}y^{\prime}.$
(11b)

Using these expressions, the no-slip boundary condition (1d) can now be
reformulated as two linear constraints on $\widehat{\omega}_{k}$:

$\forall
k,B^{+}_{k}\left(\widehat{\omega}_{k}\right)=B^{-}_{k}\left(\widehat{\omega}_{k}\right)=0,$
(12)

where

$B^{+}_{k}(f)=\int_{0}^{1}e^{-|\kappa|y}f(y)\mathrm{d}y,\quad
B^{-}_{k}(f)=\int_{0}^{1}e^{-|\kappa|(1-y)}f(y)\mathrm{d}y=0,$ (13a) for
$k\neq 0$, and
$B^{+}_{0}(f)=\int_{0}^{1}yf(y)\mathrm{d}y=0,\quad\mathrm{and}\quad
B^{-}_{0}(f)=\int_{0}^{1}(1-y)f(y)\mathrm{d}y=0.$ (13b)

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 (3), which leads to

$B^{\pm}_{k}(\partial_{y}^{2}\widehat{\omega}_{k})=\nu^{-1}B^{\pm}_{k}\left(\widehat{\boldsymbol{\nabla}(\mathbf{u}\omega)}_{k}\right).$
(14)

The above analysis ensures that the system of equations (3), (4), (5), (11)
and either one of (12) or (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 (1904).
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
$\mathrm{Re}\to\infty$:

$\omega(x,y)=\omega_{E}(x,y)+\nu^{-1/2}\omega_{P}(x,\nu^{-1/2}y)-\nu^{-1/2}\omega_{P}(x,\nu^{-1/2}(1-y))+\omega_{R}(x,y),$
(15)

where $\omega_{E}(x,y)$ is a smooth function on $\Omega\times]0,T[$,
$\omega_{P}(x,y_{P})$ is a smooth function on
$C=\mathbb{T}\times]0,\infty[\times]0,T[$ which decays rapidly when
$y_{P}\to\infty$. The indices ${}_{E},_{P}$ and R denote respectively the
Euler, Prandtl and remainder terms, and $y_{P}=\nu^{-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 $\omega_{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
$\omega_{E}$ should satisfy the incompressible Euler equations in $\Omega$,
and $\omega_{P}$ the Prandtl equations,

$\displaystyle\partial_{t}\omega_{P}+\boldsymbol{\nabla}(\mathbf{u}_{P}\omega_{P})=\partial_{y_{P}}^{2}\omega_{P}$
(16a) $\displaystyle\omega_{P}(x,y_{P},0)=0$ (16b)
$\displaystyle\psi_{P}(x,y_{P},t)=\int_{0}^{y_{P}}\mathrm{d}y_{P}^{\prime}\int_{0}^{y_{P}^{\prime}}\mathrm{d}y_{P}^{\prime\prime}\omega_{P}(x,y_{P}^{\prime\prime},t)$
(16c)
$\displaystyle\partial_{y_{P}}\omega_{P}(x,0,t)=-\partial_{x}p_{E}(x,0,t),$
(16d)

where $p_{E}$ is the pressure field calculated from $\omega_{E}$. It is
instructive to rederive the classical Neumann condition (16d) as follows.
First, by replacing $\omega$ according to (15) in (12), one obtains:

$\int_{0}^{1}\widehat{\omega}_{Ek}(y)\exp(-|\kappa|y)\mathrm{d}y+\nu^{-1/2}\int_{0}^{1}\widehat{\omega}_{Pk}(y)\exp(-|\kappa|y)\mathrm{d}y\,=\,0\,,$
(17)

and by expressing the second integral with respect to $y_{P}$ and keeping the
lowest order term in $\nu$:

$\int_{0}^{\infty}\widehat{\omega}_{Pk}(y_{P})\mathrm{d}y_{P}=-\int_{0}^{1}\widehat{\omega}_{Ek}(y)\exp(-|\kappa|y)\mathrm{d}y$
(18)

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

$-\partial_{t}\int_{0}^{1}\widehat{\omega}_{Ek}(y)\exp(-|\kappa|y)\mathrm{d}y=\int_{0}^{\infty}\partial_{y_{P}}^{2}\widehat{\omega}_{Pk}(y)=-\partial_{y_{P}}\widehat{\omega}_{Pk}(0)$
(19)

where, from the considerations in the preceding paragraph, it appears that the
left-hand side can be identified with the pressure gradient
$\partial_{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,
nothing prevents in general that the vorticity gradient in the $x$ direction
grows indefinitely, hence the possibility of finite time singularity. More
precisely, the mechanism proposed by Van Dommelen & Cowley (1990) 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^{*}$
the time at which this first occurs, and by $x^{*}$ the corresponding $x$
location. From the scaling exponents computed by Van Dommelen & Cowley (1990),
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

$k_{C}\propto(t^{*}-t)^{-3/2}.$ (20)

###  2.3 Interactive boundary layer model

As the singularity builds up in the Prandtl solution, the corresponding
Navier–Stokes solution adopts a quite different behavior. As first explained
by Elliott et al. (1983), 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

$t^{*}-t=O\left(\mathrm{Re}^{-2/11}\right).$ (21)

A rigorous asymptotic description of this new effect would require the
modification of the vorticity ansatz (15) with new BLs, both in $x$ and in
$y$, coming into play. To avoid such complications, we follow Peridier et al.
(1991b) and consider the finite Reynolds number description called 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 (15) therefore remains valid, except
that $\omega_{P}$ is replaced by $\omega_{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 (12) on $\omega$. Starting again from (17), but expanding
the exponential up to order $\mathrm{Re}^{-1/2}$, yields

$\int_{0}^{1}\widehat{\omega}_{Ik}(y_{P})\mathrm{d}y_{P}-\nu^{1/2}|\kappa|\int_{0}^{\infty}y_{P}\widehat{\omega}_{Ik}(y_{P})\mathrm{d}y_{P}=-\int_{0}^{1}\widehat{\omega}_{Ek}(y)\exp(-|\kappa|y)\mathrm{d}y$
(22)

and, following the same procedure as above, leads to a perturbed boundary
condition for $\omega_{I}$:

$\partial_{y_{P}}\widehat{\omega}_{Ik}(0)=-\iota\kappa\widehat{p}_{Ek}(0)-\nu^{1/2}|\kappa|\partial_{t}\widehat{\beta}_{Ik},$
(23)

where

${\beta}_{Ik}(x)=\int_{0}^{\infty}y_{P}{\omega}_{I}(x,y_{P})\mathrm{d}y_{P}.$
(24)

On the other hand, by multiplying (16a) by $y_{P}$ and integrating over
$y_{P}$, we obtain an expression for $\partial_{t}\beta_{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 occurence 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 (1929) to obtain an elegant asymptotic
description of the modes now known as Tollmien–Schlichting waves, and of the
corresponding stability region in the $(\mathrm{Re},k)$ plane, which was later
confirmed experimentally by Schubauer & Skramstad (1947). For a more recent
review on the subject, see Reed et al. (1996). Although most of the material
presented here is classical (see Lin (1967)), 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
$\mathrm{Re}\to\infty$ limit, which concerns us here.

For small perturbations $\delta\psi(x,y_{P},t_{2})=\phi(y_{P})e^{i(\kappa
x-\alpha t_{2})}$ to the stream function, the profile function $\phi$
satisfies the Orr–Sommerfeld equation

$(u_{P}-c)(\phi^{\prime\prime}-\nu\kappa^{2}\phi)-u_{P}^{\prime\prime}\phi=\frac{1}{i\kappa}(\phi^{\prime\prime\prime\prime}-\nu\kappa^{2}\phi^{\prime\prime}+\nu^{2}\kappa^{4}\phi),$
(25)

where $c=\alpha/\kappa$ 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

$L^{-1}\ll k\ll L^{-1}\mathrm{Re}^{1/2},$ (26)

(25) simplifies to

$(u_{P}-c)\phi^{\prime\prime}-u_{P}^{\prime\prime}\phi=\frac{1}{i\kappa}\phi^{\prime\prime\prime\prime}.$
(27)

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

$\phi(y_{P})\underset{y_{P}\to\infty}{=}A(1-\kappa\nu^{1/2}y_{P})+o(\kappa\nu^{1/2}y_{P}).$
(28)

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

####  2.4.2 Inviscid mode

Neglecting the viscous contribution in (27), we obtain

$(u_{P}-c)\phi^{\prime\prime}-u_{P}^{\prime\prime}\phi=0$ (29)

which admits the obvious regular solution

$\phi_{1,c}=u_{P}-c,$ (30)

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

$(u_{P}-c)f^{\prime\prime}+2u_{P}^{\prime}f^{\prime}=0$ (31)

and thus

$\phi_{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)\,,$
(32)

where $u_{\infty}$ is the velocity outside the boundary layer. By combining
$\phi_{1,c}$ and $\phi_{2,c}$, a solution $\phi_{out}$ satisfying the
condition (28) at $+\infty$ is readily obtained:

$\phi_{out}=\phi_{1,c}-\kappa\nu^{1/2}\phi_{2,c}\,.$ (33)

####  2.4.3 Viscous correction

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

$u_{P}^{\prime}(z_{0}(c))(y_{P}-z_{0}(c))\phi^{\prime\prime}-u_{P}^{\prime\prime}(0)\phi(0)=\frac{1}{i\kappa}\phi^{\prime\prime\prime\prime},$
(34)

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
$\eta=\varepsilon\left(y_{P}-z(c)\right)\left|\kappa
u_{P}^{\prime}(z(c))\right|5^{1/3}$, where
$\varepsilon=\mathrm{sign}(u_{P}^{\prime}(z(c)))$, leading to with $\kappa\gg
1$ to

$\eta\widetilde{\phi}^{\prime\prime}=-i\widetilde{\phi}^{\prime\prime\prime\prime},$
(35)

We are interested in a solution of this equation which remains bounded and
whose derivative tends to zero when $y_{P}\to\infty$. When $\varepsilon>0$,
this limit is equivalent to $\eta\to\infty$, and a solution to the problem was
given by Tollmien (1929) in terms of Hankel functions:

$\widetilde{\phi}_{in}(\eta)=\int_{\eta}^{\infty}\mathrm{d}\eta^{\prime}\int_{\eta^{\prime}}^{\infty}\mathrm{d}\eta^{\prime\prime}{\eta^{\prime\prime}}^{1/2}H^{(1)}_{1/3}\left(\frac{2}{3}(i\eta^{\prime\prime})^{3/2}\right)$
(36)

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

$\phi_{in}(y_{P})=\widetilde{\phi}_{in}\left(\left(y_{P}-\frac{c}{u_{P}^{\prime}(0)}\right)\left|\kappa
u_{P}^{\prime}(0)\right|^{1/3}\right)$ (37)

In the case $\varepsilon<0$, $y_{P}\to\infty$ corresponds to $\eta\to-\infty$,
and the solution $\widetilde{\phi}_{in}(\eta)$ should be adapted accordingly.

####  2.4.4 Construction of unstable modes

Now in the asymptotic regime, any admissible solution $\phi$ of (27) can
approximately be expressed as

$\phi=A\phi_{out}+B\phi_{in},$ (38)

so that the boundary conditions $\phi(0)=\phi^{\prime}(0)=0$ translate to the
linear system

$\begin{cases}A\phi_{out}(0)+B\phi_{in}(0)=0\\\
A\phi^{\prime}_{out}(0)+B\phi^{\prime}_{in}(0)=0\end{cases}.$ (39)

In order for a nontrivial solution to exist, this system should be
degenerated, i.e.

$\frac{\phi_{in}(0)}{\phi_{in}^{\prime}(0)}=\frac{\phi_{out}(0)}{\phi_{out}^{\prime}(0)}.$
(40)

On the one hand, denoting

$\eta_{w}=-c\kappa^{1/3}\left|u_{P}^{\prime}(0)\right|^{-2/3}$ (41)

the position of the wall in terms of the $\eta$ variable, it is shown
following Tollmien (1929) that

$\frac{\phi_{in}(0)}{\phi_{in}^{\prime}(0)}=-\frac{c}{u_{P}^{\prime}(0)}F(-\eta_{w})\mathrm{\quad
if\,}u_{P}^{\prime}(0)>0,$ (42)

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 (37) using quadrature formulas.

Figure 1:  Imaginary vs. real part of the Tietjens function $F(\eta)$ for
$\eta\in[2,20]$. The sampling points (every $0.1$) are indicated by crosses.

A graphical representation is given in Fig. 1.

On the other hand, denoting for convenience of notation

$I_{c}=-c\int_{0}^{\infty}\left(\frac{1}{(u_{P}(y)-c)^{2}}-\frac{1}{(u_{\infty}-c)^{2}}\right)\mathrm{d}y$
(43)

and using (33), it is shown that

$\frac{\phi_{out}(0)}{\phi_{out}^{\prime}(0)}=-\frac{c}{u_{P}^{\prime}(0)}E(\kappa\nu^{1/2},c),$
(44)

where

$E(\kappa\nu^{1/2},c)=1-\frac{\kappa\nu^{1/2}(u_{\infty}-c)^{2}}{cu_{P}^{\prime}(0)\left(1-\frac{\kappa\nu^{1/2}(u_{\infty}-c)^{2}}{c}(I_{c}-\frac{1}{u^{\prime}(0)})\right)}.$
(45)

Injecting (42) and (44) back into the degeneracy condition (40), we obtain the
dispersion relation

$F(-\eta_{w})=E(\kappa\nu^{1/2},c).$ (46)

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

The asymptotic behavior 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

$I_{c}\underset{c\to
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)-ci\pi\Delta_{u}+cK_{R}+o(c)$
(47)

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

$\Delta_{u}=\sum_{\\{y_{P}>0\mid
u_{P}(y_{P})=0\\}}\frac{u_{P}^{\prime\prime}(y_{P})}{u_{P}^{\prime}(y_{P})^{3}}.$
(48)

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

$\Delta^{\pm}_{u}=\begin{cases}\Delta_{u}\mathrm{\,\,if\,\,}\Re(c/u^{\prime}_{P}(0))<0\\\
\pm\frac{u_{P}^{\prime\prime}(0)}{u_{P}^{\prime}(0)^{3}}+\Delta_{u}\mathrm{\,\,if\,\,}\Re(c/u^{\prime}_{P}(0))>0\end{cases}.$
(49)

By injecting our estimate for $I_{c}$ in (45), and assuming for simplicity
that $\kappa\nu^{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:

$\displaystyle\Re(E)$ $\displaystyle\sim
1-\frac{\kappa\nu^{1/2}u_{\infty}^{2}}{cu_{P}^{\prime}(0)}$ (50a)
$\displaystyle\Im(E)$ $\displaystyle\sim\frac{\kappa^{2}\nu
u_{\infty}^{4}\Delta^{\pm}_{u}\pi}{cu_{P}^{\prime}(0)},$ (50b)

where we have assumed that $\Delta^{\pm}_{u}\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, (40) can be satisfied
only in the neighborhood of points where $F$ is purely real. This happens for
$\eta_{w}\to-\infty$, as well as for a certain finite value
$\eta_{0}\simeq-2.3$ where $F(-\eta_{0})\simeq 0.56$ (see Fig. 1).

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

$k_{1}\simeq 0.161\,|u_{P}^{\prime}(0)|^{5/4}u_{\infty}^{-3/2}\nu^{-3/8},$
(51)

beyond which the modes are unstable. In the case $u^{\prime}_{P}(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 $\eta_{w}\to-\infty$, corresponds to $F(-\eta_{w})\to
0$, so from (50a)

$\frac{\kappa\nu^{1/2}u_{\infty}^{2}}{cu_{P}^{\prime}(0)}\simeq 1.$ (52)

To obtain another equation relating $c$ and $\kappa$, we use the estimate
$F(-\eta_{w})\underset{\eta\to-\infty}{\sim}-e^{i\pi/4}|\eta|^{-3/2},$ which,
combined with (50b), leads to

$\frac{\kappa^{2}\nu
u_{\infty}^{4}\Delta^{\pm}_{u}\pi}{cu_{P}^{\prime}(0)}\simeq\frac{\sqrt{2}}{2}|\eta_{w}|^{-3/2}.$
(53)

This equation has a simple root iff $\Delta_{u}>0$ and $u^{\prime}_{P}(0)<0$,
in which case we obtain a second critical wavenumber

$k_{2}=0.0968\,|u_{P}^{\prime}(0)|^{1/2}|\Delta_{u}|^{-1/3}|u_{\infty}|^{-5/3}\nu^{-5/12}$
(54)

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\nu^{-1/2}$, which
implies that it extends at least up to wavenumbers scaling like $\nu^{-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 further down.

Another important quantity we need to estimate is the growth rate of the
unstable modes. From (45), we see that since $c$ remains small in absolute
value, the growth rate $\alpha$ of the instability satisfies

$\alpha\sim\frac{u_{\infty}^{2}\Im(F(-\eta_{w}))}{|u^{\prime}_{P}(0)|}\nu^{1/2}\kappa^{2}$
(55)

To sum up, the generic instability expected to play a role in such boundary
layer flows in the inviscid limit manifest itself by the growth of wave
packets in the vicinity of the boundary confined in physical space to regions
where $u^{\prime}_{P}(0)<0$ (recirculation bubbles), and whose parallel
wavenumber support extends from $O(\mathrm{Re}^{3/8})$ to at least
$O(\mathrm{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 (1965). Going
back to the general expression (33) for the outer solution, we obtain in the
case $u_{P}^{\prime}(0)=0$ that

$E(\kappa\nu^{1/2},c)=\frac{\phi_{out}(0)}{\phi_{out}^{\prime}(0)}=-\frac{c^{2}}{\kappa\nu^{1/2}(u_{\infty}-c)^{2}}-cI_{c}$
(56)

or, with (87), and when $c$ is close to the real axis:

$\displaystyle\Re(E)$
$\displaystyle\simeq-\frac{c^{2}}{\kappa\nu^{1/2}u_{\infty}^{2}}$ (57a)
$\displaystyle\Im(E)$
$\displaystyle\simeq-\frac{\pi}{4}\left(\frac{2c}{u_{P}^{\prime\prime}(0)}\right)^{1/2}$
(57b)

Concerning the inner solution, (42) is replaced by

$\frac{\phi_{in}(0)}{\phi^{\prime}_{in}(0)}=-z(c)F(-\eta_{w}),$ (58)

which, combined with (57b), yields

$\Im(F(-\eta_{w}))=-\frac{\pi}{4},$ (59)

or equivalently, according to Hughes & Reid (1965),

$\eta_{w}\simeq-0.488,\quad\Re(F(-\eta_{w}))\simeq 1.580,$

and therefore

$\displaystyle\frac{c^{3/2}}{\kappa\nu^{1/2}u_{\infty}^{2}}=\sqrt{\frac{2}{u_{P}^{\prime\prime}(0)}}1.580,$
(60)
$\displaystyle\frac{2c}{u_{P}^{\prime\prime}(0)}\kappa^{1/3}(2cu^{\prime\prime}(0))^{1/3}=0.488.$
(61)

which finally gives us the critical wavenumber

$k_{1}=0.0279\,u_{P}^{\prime\prime}(0)^{10/11}u_{\infty}^{16/11}\nu^{-4/11}$

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

$\alpha\sim u_{\infty}\Im(F(-\eta_{w}))\nu^{1/4}\kappa^{3/2}$ (62)

####  2.4.6 Physical interpretation

In this section we formulate some conjectures relevant to 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 $\mathrm{Re}\to\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.
Even though the range $[k_{1},k_{2}]$ is already subject to an instability,
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 is 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 $\kappa$ by the characteristic excitation (20) generated by
the Prandtl dynamics some time $t^{*}-t$ before the singularity, we obtain

$\alpha\propto\nu^{1/2}(t^{*}-t)^{-3}.$ (63)

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

$t^{*}-t=O(\mathrm{Re}^{-1/4}).$ (64)

By comparing this result with (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 of 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. (2014), who pointed out
that what they call a large scale interaction always preceeds the approach to
detachment.

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

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

##  3 Solvers

###  3.1 Setup

To trigger an unsteady separation process, we have chosen an initial
configuration inspired by the dipole of Orlandi (1990), later modified by
Clercx & van Heijst (2002). 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 $3/4$ of the domain size for a given
$\mathrm{Re}$. It is defined in terms of its stream function as follows:

$\displaystyle\psi_{i}(x,y)$
$\displaystyle=Axy\exp\left(-\frac{(x-x_{0})^{2}+(y-y_{0})^{2}}{2s^{2}}\right),$
(65a)

where $A=0.625847306637464$ determines amplitude of the vortices,
$s=6.3661977236758$ 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

$v_{i}(x,y=0)\approx 10^{-15}.$

which is anyway of the same order as the round-off error in double precision
arithmetics.

Due 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 $\mathbf{u}_{i}$ in $K$ are shown in Fig. 2. The definitions
and initial values of several integral quantities which we will be important
in our study are given in Table 1.

quantity | enstrophy | maximum vorticity | energy | maximum velocity
---|---|---|---|---
definition | $\displaystyle{\frac{1}{2}\int_{\Omega}\omega^{2}}$ | $\displaystyle{\max_{\Omega}\omega}$ | $\displaystyle{\frac{1}{2}\int_{\Omega}|\mathbf{u}|^{2}}$ | $\displaystyle{\max_{\Omega}|\mathbf{u}|}$
initial value | $7.48\cdot 10^{-3}$ | $1.00$ | $9.05\cdot 10^{-6}$ | $2.42\cdot 10^{-2}$
Table 1:  Definitions and initial values of several integral quantities of
interest. Figure 2:  Streamlines of initial velocity field in the subdomain
$K=[1/2,1]\times[0,1/2]$.

In this work we shall analyze the flows obtained by solving the Navier–Stokes
equations numerically up to $T=57.05$ for $\nu$ decreasing from $4\cdot
10^{-7}$ to $5\cdot 10^{-8}$ by factors of $\sqrt{2}$ (i.e., $7$ different
values in total). Reynolds numbers corresponding to these values of $\nu$ are
defined according to (2), where $U\simeq 2.42\cdot 10^{-2}$ is the initial
maximum velocity, and $L=2s\simeq 1.27\cdot 10^{-1}$ is a measure of the size
of the quadrupole. Both $\nu$ and $\mathrm{Re}$ are provided up to 3
significant digits in Table 2. To facilitate comparison with previous results
concerning dipole-wall collisions, we have also included the Reynolds number
$\mathrm{Re}_{rms}$ computed from the RMS velocity $U_{rms}=4.25\cdot 10^{-3}$
and channel width instead, which is of the same order of magnitude.

case | I | II | III | IV | V | VI | VII | VIII | IX
---|---|---|---|---|---|---|---|---|---
$10^{8}\nu$ | $40.0$ | $28.3$ | $20.0$ | $14.1$ | $10.0$ | $7.07$ | $5$ | $3.54$ | $2.5$
$10^{-2}\mathrm{Re}_{rms}$ | $106$ | $150$ | $213$ | $301$ | $425$ | $601$ | $851$ | $1200$ | $1700$
$10^{-2}\mathrm{Re}$ | $76.8$ | $109$ | $154$ | $217$ | $308$ | $435$ | $615$ | $870$ | $1230$
$10^{3}\Delta_{y_{E},1}$ | $1.90$ | $1.91$ | $1.92$ | $1.92$ | $1.93$ | $1.93$ | $1.93$ | $1.94$ | $1.94$
$10^{5}\Delta_{y_{P},1}$ | $15.4$ | $13.0$ | $11.0$ | $9.25$ | $7.80$ | $6.57$ | $5.54$ | $4.66$ | $3.93$
$10^{4}\Delta_{x,1}$ | $9.77$ | $9.77$ | $9.77$ | $9.77$ | $9.77$ | $4.88$ | $4.88$ | $4.88$ | $4.88$
$10^{3}\Delta_{y_{E},2}$ | $1.76$ | $1.87$ | $1.95$ | $2.02$ | $2.07$ | $2.12$ | $2.16$ | $2.19$ | $2.22$
$10^{5}\Delta_{y_{P},2}$ | $14.4$ | $12.7$ | $11.2$ | $9.71$ | $8.39$ | $7.22$ | $6.18$ | $5.28$ | $4.50$
$10^{5}\Delta_{y_{K},2}$ | $11.4$ | $8.46$ | $6.24$ | $4.56$ | $3.32$ | $2.40$ | $1.73$ | $1.24$ | $0.889$
$10^{4}\Delta_{x,2}$ | $9.77$ | $9.77$ | $4.88$ | $2.44$ | $1.22$ | $0.610$ | $0.305$ | $0.305$ | $0.153$
Table 2:  Parameters of numerical experiments. All figures in this table are
given up to 3 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. Since the $y_{P}$
direction is not periodic, derivatives in the $y_{P}$ direction have to be
treated differently. The Chebychev scheme is accurate but very costly, and
also imposes Gauss collocation points which are not optimal for our problem.
We have therefore prefered to turn to fifth order compact finite differences
(Lele, 1992; Gamet et al., 1999). Denoting by $(f_{i})$ approximate values of
a function $f$ on the uniform grid defined by
$y_{1,i}=\frac{i}{N_{y_{P}}-1}(0\leq i<N_{y_{P}})$, and $(f_{i}^{\prime})$,
$(f_{i}^{\prime\prime})$ approximations of its first and second derivatives at
the same locations, we impose fifth-order accuracy by requesting that

$\displaystyle\alpha_{i}f_{i-1}^{\prime}+f_{i}^{\prime}+\beta_{i}f_{i+1}^{\prime}$
$\displaystyle=A_{i}f_{i-1}+B_{i}f_{i}+C_{i}f_{i+1}$ (66a)
$\displaystyle\gamma_{i}f_{i-1}^{\prime\prime}+f_{i}^{\prime\prime}+\delta_{i}f_{i+1}^{\prime\prime}$
$\displaystyle=D_{i}f_{i-1}+E_{i}f_{i}+F_{i}f_{i+1}$ (66b)

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\leq 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
$\partial_{y}(v\omega)$, 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 $\partial_{y}^{2}\omega$, they should follow from the
boundary conditions (14). To derive them, the integrals are first discretized
by a fifth-order local quadrature formula. To preserve accuracy, $\omega$ is
expanded locally into its Taylor polynomial form, and the contribution of the
$k$-depdendent 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 (Li et al., 1999; Demmel et al., 1999), and the PetscC
library (Balay et al., 2013) is used for matrix arithmetics. Note that due to
the dependency of the integral constraint on $k$, the number of LU
factorizations is multiplied by $K$. The cost of these factorizations is
considerable, and they are tractable only under the condition that $N_{y}$ is
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 $\Delta_{y_{P},1}$ between 0 and

$L_{1}=166\,\mathrm{Re}^{-1/2}L,$

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 $\Delta_{y_{E},1}$.

Figure 3:  Schematics of discretization grids used by the Navier–Stokes solver
for $t\leq 54$ (left) and $t>54$ (right). 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 $\mathrm{Re}$. The new grid has $N_{y,2}=449$ points in the
$y$ direction. The grid spacing is set to $\Delta_{y_{K},2}$ between 0 and

$L_{2}=1210\,\mathrm{Re}^{-1}L,$

then to $\Delta_{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
Fig. 3.

###  3.3 Prandtl–Euler solver

Following previous work (see e.g. Nguyen van yen et al., 2009), the Euler
equation is approximated by the Navier–Stokes equations with hyperdissipation,
i.e., the dissipation term $\nu\Delta\omega$ is replaced by
$-\nu_{2}(-\Delta)^{2}\omega$. This approximation is second order in space
(Kato, 1972), 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=4096$ grid points on the whole domain. A
low-storage third order Runge–Kutta scheme is employed for time
discretization, the time step being adjusted dynamically to satisfy the CFL
condition. The hyperviscosity parameter $\nu_{2}$ was set to $2.146\cdot
10^{-13}$, which was found to sufficiently regularize the solution.

To solve the initial value problem for the Prandtl equations (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
solution does not depend on its value on the time interval considered. 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}\partial_{y_{P}}\omega_{P}$, the
equations at the edges are obtained by shifting the stencils so that they
remain inside the computational grid (no additional condition is included).
For the dissipation term $\partial_{y_{P}}^{2}\omega_{P}$, the integral
constraint (19) is rewritten as follows:

$\int_{0}^{\infty}\partial_{y_{P}}^{2}\omega_{P}(x,y_{P},t)\mathrm{d}y_{P}=\partial_{x}p_{E}(x,0,t),$
(67)

which is enforced as for the Navier–Stokes solver. Finally, the system is
closed by imposing that $\omega_{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 (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

$\Delta t<C\frac{\Delta_{x}^{2}}{\sqrt{\nu}},$

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 indepentently 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 $\omega(L_{y_{P}})=0$ and
$\phi^{\prime}(L_{y_{P}})+k\nu^{1/2}\phi(L_{y_{P}})=0$, which follow from the
reconnection with a potential solution at large $y_{P}$.

A second order finite differences 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 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 behavior of the various solutions well before the Prandtl singularity time
is well understood and we present it only for the sake of completeness.

Figure 4:  Top: contour lines of vorticity field at $t=30$ (left) and $t=50$
(right), for $\omega=0.1,0.2,\ldots,1$. Bottom: contour lines of corresponding
boundary layer vorticity fields in Prandtl variables, for
$\omega=-|\omega|_{max}/6,-2|\omega|_{max}/6,\ldots,-5|\omega|_{max}/6$, where
$|\omega|_{max}=3.75\cdot 10^{-4}$ and $|\omega|_{max}=8.70\cdot 10^{-3}$,
respectively.

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 Fig. 4 (top, left). The Navier–Stokes flow in the
Prandtl boundary layer units (bottom, left) 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 hand
(Fig. 4, right). However, as expected, the discrepancy between Prandtl and
Navier–Stokes flows has increased.

Figure 5:  Contour lines of vorticity field in Prandtl variables at $t=54$
(left) and $t=55.3$ (right), for
$\omega=5|\omega|_{max}/6,\ldots,-|\omega|_{max}/6,0,-|\omega|_{max}/6,-2|\omega|_{max}/6,\ldots,-5|\omega|_{max}/6$,
where $|\omega|_{max}=1.25\cdot 10^{-2}$ and $|\omega|_{max}=1.30\cdot
10^{-2}$, respectively. The $\omega=0$ contour line is shown in bold.

At $t=54$ (Fig. 5, left) 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

First signs of a qualitatively different behavior become visible shortly
thereafter, as shown for example at $t=55.3$ in Fig. 5 (right). The contour
lines of the Prandtl vorticity have become very concentrated around
$x^{*}=0.556$, indicating the formation of a finite time singularity with
precisely the qualitative features predicted by Van Dommelen & Cowley (1990),
in particular a blow-up of the wall-normal velocity associated to an infinite
acceleration of fluid particles away from the wall.

As the Prandtl solution approaches its singularity time $t^{*}\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. (1983), the Navier–Stokes solution departs from
the Prandtl behavior when the potential flow perturbation due to the presence
of the boundary layer starts to perturb the wall pressure gradient.

(a) $t=50$ (b) $t=54$ (c) $t=50$ (d) $t=54$

Figure 6:  Comparison between Navier–Stokes (black), Prandtl–Euler (blue) and
interactive boundary layer (red) models for varying Reynolds numbers, at
$t=50$ (left) and $t=54$ (right).

By comparing the Navier–Stokes, Prandtl and interactive boundary layer models
at different Reynolds numbers for $t=50$, we observe good agreement between
all models (Fig. 6 (a)). For $t=54$ (Fig. 6 (b)) it can be noted that the IBL
solution has indeed slightly departed away from the Prandtl solution, but to a
point which fails way short of capturing the full behavior of the NS solution.
This effect, which corresponds in principle to the large scale interaction
also described by Gargano et al. (2014), 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$.

(a) $t=54$ (b) $t=55.3$

Figure 7:  Boundary vorticity for varying Reynolds number.

(a) $t=54$ (b) $t=55.3$

Figure 8:  Fourier spectra of boundary vorticity for varying Reynolds number.
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 $\mathrm{Re}$ at the instants
$t=54$ and $t=55.3$ in Fig. (7), we observe that for increasing $\mathrm{Re}$
the elbow structure looks more and more like a wave-packet 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 Fig. 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 $\omega$
along the wall. The NSE solutions spectra all develop a distinctive bump in a
wavenumber range. Both the width and $k$ location of this bump increase with
Reynolds and with time. Interestingly, for the largest $\mathrm{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 around $k=20$. This confirms a posteriori the validity of the slow-
varying flow approximation used in deriving the asymptotic stability results
of Sec. 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=123075$. The upper and lower
marginal wavenumbers estimated from our asymptotic analysis (Eq. 51 and 54)
are shown in dashed lines.

We now would 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 (Fig. 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 Fig. 6 (b). In $k$-space, the Orr–Sommerfeld computations predict that
the instability should start around $k=30$ for the high $\mathrm{Re}$
considered, which is in very good agreement as well with the wavenumber at
which the corresponding spectrum in Fig. 8 (a) starts to exceed the reference
Prandtl solution (shown in dashed lines). This effect becomes more pronounced
at $t=55.3,$ shown in Fig. 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^{\prime}_{P}(0)=0$ is approached. Nevertheless, the overall instability
region is qualitatively well captured by the criterion $u^{\prime}_{P}(0)<0$.

###  4.4 Detachment and production of dissipative structures

Figure 10:  Contour lines of vorticity field in Prandtl variables at $t=56$
(left) and $t=56.9$ (right).

The instability process which we have seen at play above introduces a new
vorticity scaling, $\mathrm{Re}^{1}$, very close to the wall. This new scaling
is difficult to notice at first, because it is hidden behind the preexisting
large negative vorticity of the boundary layer.

Figure 11:  Reynolds number scalings of maximum vorticity (left) and enstrophy
(right).

A simple trick to observe it more easily is to consider only the maximum
vorticity of positive sign (Fig. 11, left). This quantity scales like
$\mathrm{Re}^{1/2}$ at $t=53$, and at $t=56.9$ it has clearly transited to the
stronger $\mathrm{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, as illustrated in Fig. 10.

###  4.5 Later evolution

Figure 12:  Contour lines of vorticity field at $t=80$.

At much later times, the Euler and Navier–Stokes solutions have become
entirely different (Fig. 12). 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 lead 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 practices in numerical fluid dynamics, we have taken care of using
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 13:  Comparison between Navier–Stokes solutions at the same Reynolds
number with two numerical resolutions, in order to check for convergence.
Left: $t=56$. Right: $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 Fig. 13, the computation used
in our analysis and the one obtained with twice the grid spacing agree well.

##  5 Discussion

There are several features of the numerical solutions which had not been
observed in previous work. The most striking one is the appearance of the
scaling $\mathrm{Re}^{1}$ for the vorticity maximum, which takes precedence,
at the singularity time, over the weaker Prandtl scaling $\mathrm{Re}^{1/2}$.
Even more strikingly, as seen on the graphs of the 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 (12) on the vorticity field, which are themselves
consequences of the no-slip boundary conditions combined with the
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 (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 $\mathrm{Re}^{1/2}$, while
the amplitude of vorticity continues to scale as $\mathrm{Re}^{1/2}$ due to
the presence of a Prandtl boundary layer. Therefore, as soon as the buildup of
the Prandtl singularity sufficiently excites those wavenumbers, their
superposition induces a $\mathrm{Re}^{1}$ scaling for the amplitude of
$\omega$. In the linear phase, the thickness of the corresponding new wall-
normal sublayer scales like $\mathrm{Re}^{2/3}$, but as soon as the
instability becomes nonlinear, vorticity transport induces excitation of
scales as fine as $\mathrm{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 modeled 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 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, 1991). 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

$\langle
U(y)\rangle\simeq\frac{U_{\tau}}{K_{\mathrm{karman}}}\log\left(\frac{yU_{\tau}}{\nu}\right),$
(68)

where

$U_{\tau}=\sqrt{\nu\left\langle\left.\frac{\mathrm{d}U}{\mathrm{d}y}\right|_{y=0}\right\rangle}$
(69)

is the so-called friction velocity. This behavior is confirmed by the most
recent experiments, with subtle corrections. An important consequence is that
the bulk velocity and $U_{\tau}$ have the same scaling with $\mathrm{Re}$ up
to a logarithmic factor. Then, from (69), one can immediately deduce that
$\left.\frac{\mathrm{d}\langle u\rangle}{\mathrm{d}y}\right|_{y=0}$ scales
like $\mathrm{Re}^{1}$ up to a logarithmic factor, which can be seen as the
statistical signature of the existence of a boundary layer of thickness
$\mathrm{Re}^{-1}$ in the neighborhood 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_{\tau}$ scales with
the bulk velocity, and this scaling was introduced by von Kármán (1921)
precisely to account for the behavior of the drag coefficient at high Reynolds
number, which was recognized as the essential issue in 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 Claude Bardos, David Gérard-Varet, Helena
Nussenzveig Lopes, Milton Lopes Filho, and especially Stephen Cowley for
important discussions. RNVY thanks the Humboldt Foundation for supporting his
research by a post-doctoral fellowship. MF and KS acknowledge support by the
French Research Federation for Fusion Studies within the framework of the
European Fusion Development Agreement (EFDA).

##  Appendix A Estimation of the integral $\boldmath I_{c}$

In this appendix, we establish the estimate (47) for the integral $I_{c}$
defined by (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
behavior of $I_{c}$ when $c\to 0$ is controlled by the behavior of $u$ around
its zeros on $[0,+\infty[$. Hence, let $(y_{i})_{1\leq i\leq 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 neighborhood of
the real axis, which ensures that there exists complex neighborhoods $Y_{i}$
of $y_{i}$ such that $u$ is a local holomorphism $u:Y_{i}\to u(Y_{i})$, and
$u(Y_{i})$ are neighborhoods 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

$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))},$
(70)

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 equiintegrable on $Y_{i}$. Therefore, if we define $\xi_{0}=0$,
$\xi_{i}=\frac{1}{2}\left(y_{i+1}+y_{i}\right)$ for $i<m$ and
$\xi_{m}=2y_{m}+1$, $f_{i}(y,c)$ is equiintegrable on $[\xi_{i},\xi_{i+1}]$,
and by the Vitali convergence theorem, the pointwise limit

$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}}$
(71)

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

$\int_{\xi_{i}}^{\xi_{i+1}}f_{i}(y,c)\mathrm{d}y\xrightarrow[c\to
0]{}\int_{\xi_{i}}^{\xi_{i+1}}f_{i}(y,0)\mathrm{d}y$ (72)

Now letting

$\displaystyle J_{i,c}=$
$\displaystyle\int_{\xi_{i}}^{\xi_{i+1}}\frac{\mathrm{d}y}{u^{\prime}(z_{i}(c))^{2}(y-z_{i}(c))^{2}}$
(73a) $\displaystyle K_{i,c}=$
$\displaystyle-\int_{\xi_{i}}^{\xi_{i+1}}\frac{u^{\prime\prime}(z_{i}(c))\mathrm{d}y}{u^{\prime}(z_{i}(c))^{3}(y-z_{i}(c))},$
(73b)

by construction

$I_{c}=-c\sum_{i=0}^{m-1}\int_{\xi_{i}}^{\xi_{i+1}}\left(f_{i}(y,c)\mathrm{d}y+J_{i,c}+K_{i,c}\right)-c\int_{\xi_{m}}^{\infty}\left(\frac{1}{(u(y)-c)^{2}}-\frac{1}{(u_{\infty}-c)^{2}}\right)$
(74)

or, by (72) and the fact that $u$ does not vanish on $[\xi_{m},+\infty]$,

$\displaystyle I_{c}\underset{c\to 0}{=}$
$\displaystyle-c\sum_{i=0}^{m-1}\left(\int_{\xi_{i}}^{\xi_{i+1}}f_{i}(y,0)\mathrm{d}y+J_{i,c}+K_{i,c}\right)-c\int_{\xi_{m}}^{\infty}\left(\frac{1}{u(y)^{2}}-\frac{1}{u_{\infty}^{2}}\right)+o(c)$
(75a) $\displaystyle\underset{c\to 0}{=}$
$\displaystyle-c\sum_{i=0}^{m-1}\left(J_{i,c}+K_{i,c}\right)+cK_{R}+o(c)$
(75b)

so that the dominant behavior of $I_{c}$ is controlled by $J_{i,c}$ and
$K_{i,c}$, which we now proceed to estimate. Since

$J_{i,c}=-\frac{1}{u^{\prime}(z_{i}(c))^{2}}\left(\frac{1}{\xi_{i+1}-z_{i}(c)}-\frac{1}{\xi_{i}-z_{i}(c)}\right),$
(76)

$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

$z_{0}(c)\underset{c\to 0}{=}\frac{c}{u^{\prime}(0)}+o(c),$ (77)

and therefore

$J_{0,c}\underset{c\to 0}{=}-\frac{1}{u^{\prime}(0)c}+K_{R}+o(1).$ (78)

Now

$\displaystyle K_{i,c}=$
$\displaystyle-\frac{u^{\prime\prime}(z_{i}(c))}{u^{\prime}(z_{i}(c))^{3}}\left(\ln(\xi_{i+1}-z_{i}(c))-\ln(\xi_{i}-z_{i}(c))\right),$
(79)

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

$\displaystyle K_{i,c}\underset{c\to 0}{=}$
$\displaystyle-\frac{u^{\prime\prime}(y_{i})}{u^{\prime}(y_{i})^{3}}\left(\ln(\xi_{i+1}-y_{i})-\ln(\xi_{i}-y_{i})\right)+o(1)$
(80a) $\displaystyle\underset{c\to 0}{=}$
$\displaystyle\iota\pi\frac{u^{\prime\prime}(y_{i})}{u^{\prime}(y_{i})^{3}}+K_{R}+o(1),$
(80b)

whereas we obtain

$K_{0,c}\underset{c\to
0}{=}-\frac{u^{\prime\prime}(0)}{u^{\prime}(0)^{3}}\left(-\ln\left(-\frac{c}{u^{\prime}(0)}\right)\right)+K_{R}+o(1)$
(81)

Finally, combining (75b), (78), (80b) and (81), we get the desired estimate of
$I_{c}$:

$I_{c}\underset{c\to
0}{=}\frac{1}{u^{\prime}(0)}-c\frac{u^{\prime\prime}(0)}{u^{\prime}(0)^{3}}\ln\left(-\frac{c}{u^{\prime}(0)}\right)-ci\pi\sum_{i=1}^{m-1}\frac{u^{\prime\prime}(y_{i})}{u^{\prime}(y_{i})^{3}}+cK_{R}+o(c)$
(82)

If $u^{\prime}(0)=0$, we cannot apply Vitali’s theorem to (70) for $i=0$
because the second and third term diverge when $c\to 0$. Hughes & Reid (1965)
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

$(u-c)^{2}=(v-\sqrt{c})(v+\sqrt{c}),$ (83)

with

$v^{\prime}(z)^{2}\xrightarrow[z\to 0]{}\frac{u^{\prime\prime}(0)}{2},\quad
v^{\prime\prime}(z)\xrightarrow[z\to
0]{}\pm\sqrt{\frac{2}{u^{\prime\prime}(0)}}u^{\prime\prime\prime}(0),$ (84)

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

$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))}.$
(85)

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^{-}$. Due to
the order of the different terms, it is sufficient to keep only one of them

$I_{c}\underset{c\to
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}},$
(86)

and after computing the residuals we finally obtain

$I_{c}\underset{c\to 0}{=}(8u^{\prime\prime}(0)c)^{-1/2}\iota\pi+o(c^{-1/2})$
(87)

##  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 test case for the Navier–Stokes solver, the setup designed by Kramer et al.
(2007) was considered. Contrary to the quadrupole setup which has been studied
in the body of the present paper, the dipole is not symmetric with respect to
the channel centerline. The full span of the channel and the walls on both
sides therefore needs to be taken into account. Two runs with
$\mathrm{Re}=650$ and $\mathrm{Re}=2500$ were performed, respectively with
$512\times 255$ and $1024\times 511$ uniformly distributed grid points.

|  | Current results | Kramer et al. (2007)
---|---|---|---
$\mathrm{Re}$ | $t$ | $x_{d}$ | $y_{d}$ | $\omega_{d}$ | $x_{d}$ | $y_{d}$ | $\omega_{d}$
| $0.0$ | $0.0977$ | $1.0000$ | $316.70$ | $0.1000$ | $1.0000$ | $316.84$
$625$ | $0.6$ | $0.1680$ | $0.1836$ | $158.51$ | $0.1656$ | $0.1827$ | $158.70$
$625$ | $1.0$ | $0.2578$ | $0.1953$ | $102.53$ | $0.2543$ | $0.1949$ | $102.64$
$2500$ | $0.6$ | $0.1797$ | $0.1016$ | $261.83$ | $0.1654$ | $0.1042$ | $261.83$
$2500$ | $1.0$ | $0.1562$ | $0.1445$ | $231.22$ | $0.2185$ | $0.1738$ | $231.40$
Table 3:  Location $(x_{d},y_{d})$ and value $(\omega_{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.

In order to make a quantitative comparison, the same procedure used by Kramer
et al. (2007) 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 misstake 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.

###  B.2 Prandtl–Euler solver

For the Prandtl solver, the classical impulsively started cylinder studied by
VDS in the Lagrangian framework is employed as a test case. It corresponds to
a constant boundary pressure gradient given by

$\partial_{x}p_{E}(x,0,t)=-\sin(x)\cos(x),$ (88)

and an initial vorticity which is a Dirac distribution

$\omega_{P}(x,y,0)=-\delta_{0}(y)\sin(x)$ (89)

For spatial discretization, $1023$ grid points are considered on the interval
$[0,\pi]$ 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

$\omega_{P}(x_{i},0,0)=-\frac{\sin(x_{i})}{\Delta
x},\quad\omega_{P}(x_{i},y,0)=0\mathrm{\ if\ }y>0.$ (90)

The results are then compared with Fig. 10 and Table II of VDS.

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

Fig. 14 shows the wall shear stress for different time instants. It is in very
good qualitative agreement with Fig. 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^{\prime\prime}_{w}$ 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$).

###  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.

Figure 15:  Unstable region for the Blasius boundary layer.

The obtained marginal stability curve (Fig. 15) is in good agreement with
Fig.16.11 provided by Schlichting (1979).

##  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 (Fig. 16 top, left) remains very close to the
Euler vorticity field (Fig. 16 top, right). Plotting the Navier-Stokes flow in
the Prandtl boundary layer units (Fig. 16 top, left) reveals that it is
smooth, and very well approximated by the solution of the Prandtl equations
(Fig. 16 top, right). The vorticity along the boundary (Fig. 18, top, left)
converges to the Prandtl values as the Reynolds number increases.

As the dipole approaches the wall (Fig. 16, 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 boundary vorticity (Fig. 18, top, right) reveals a larger difference
between Prandtl and Navier–Stokes flows than at $t=30$.

As the Prandtl solution approaches its singularity time $t^{*}\approx 55.6$
(Fig. 16, 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 (Fig. 18, bottom, left).

After the singularity, at $t=57$ (Fig. 17, top) oscillations in the vorticity
appear, while the bulk flow still looks similar for Navier–Stokes and Euler.
Following the new vorticity extremum which 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 (Fig. 18,
bottom, right).

At much later times, $t=100$ (Fig. 17, 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 lead to the formation of two new vortices (shown in
cyan and magenta in the Fig. 17, 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.

Figure 16:  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. Figure 17:
Continuation of Fig. 16 at time instants $t=57$ and $100$. The Prandtl
solution is singular and thus not shown.

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

## References

  * Balay et al. (2013) Balay, Satish, Adams, Mark F., Brown, Jed, Brune, Peter, Buschelman, Kris, Eijkhout, Victor, Gropp, William D., Kaushik, Dinesh, Knepley, Matthew G., McInnes, Lois Curfman, Rupp, Karl, Smith, Barry F. & Zhang, Hong 2013 PETSc users manual. Tech. Rep. ANL-95/11 - Revision 3.4. Argonne National Laboratory.
  * Blasius (1908) Blasius, H. 1908 Grenzschichten in Flüssigkeiten mit kleiner Reibung. Z. Math. Phys. 56, 1–37, engl. trans. in NACA TM 1256 (1950).
  * Bowles (2006) Bowles, R. 2006 Lighthill and the triple-deck, separation and transition. Journal of Engineering Mathematics 56 (4), 445–460.
  * Bowles et al. (2003) Bowles, Robert I, Davies, Christopher & Smith, Frank T 2003 On the spiking stages in deep transition and unsteady separation. Journal of Engineering Mathematics 45 (3), 227–245.
  * Brinckman & Walker (2001) Brinckman, KW & Walker, JDA 2001 Instability in a viscous flow driven by streamwise vortices. Journal of Fluid Mechanics 432, 127–166.
  * Cassel (2000) Cassel, Kevin W. 2000 A comparison of Navier–Stokes solutions with the theoretical description of unsteady separation. Philosophical Transactions: Mathematical, Physical and Engineering Sciences 358 (1777), 3207–3227.
  * Cassel & Obabko (2010) Cassel, Kevin W. & Obabko, Aleksandr V. 2010 A Rayleigh instability in a vortex-induced unsteady boundary layer. Physica Scripta T142, 1–14.
  * Clercx & van Heijst (2002) Clercx, H. & van Heijst, G. J. 2002 Dissipation of kinetic energy in two-dimensional bounded flows. Phys. Rev. E 65, 066305.
  * Cowley (2002) Cowley, S.J. 2002 Laminar boundary-layer theory: A 20th century paradox? Mechanics for a New Mellennium, 389–412.
  * Demmel et al. (1999) Demmel, James W., Gilbert, John R. & Li, Xiaoye S. 1999 An asynchronous parallel supernodal algorithm for sparse Gaussian elimination. SIAM J. Matrix Analysis and Applications 20 (4), 915–952.
  * E & Engquist (1997) E, Weinan & Engquist, B. 1997 Blowup of solutions of the unsteady Prandtl’s equation. Communications on Pure and Applied Mathematics 50 (12), 1287–1293.
  * E & Liu (1996) E, Weinan & Liu, J.G. 1996 Vorticity boundary condition and related issues for finite difference schemes. Journal of Computational Physics 124 (2), 368–382.
  * Elliott et al. (1983) Elliott, JW, Smith, FT & Cowley, SJ 1983 Breakdown of boundary layers:(i) on moving surfaces;(ii) in semi-similar unsteady flow;(iii) in fully unsteady flow. Geophysical & Astrophysical Fluid Dynamics 25 (1-2), 77–138.
  * Gamet et al. (1999) 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. International Journal for Numerical Methods in Fluids 29 (2), 159–191.
  * Gargano et al. (2009) Gargano, F., Sammartino, M. & Sciacca, V. 2009 Singularity formation for Prandtl’s equations. Physica D 238 (19), 1975–1991.
  * Gargano et al. (2011) Gargano, F., Sammartino, M. & Sciacca, V. 2011 High Reynolds number Navier–Stokes solutions and boundary layer separation induced by a rectilinear vortex. Computers & Fluids 52, 73–91.
  * Gargano et al. (2014) Gargano, F., Sammartino, M., Sciacca, V. & Cassel, K.W., 2014. Analysis of complex singularities in high-Reynolds-number Navier–Stokes solutions. Journal of Fluid Mechanics, 747, 381–421.
  * Gérard-Varet & Dormy (2010) Gérard-Varet, D. & Dormy, E. 2010 On the ill-posedness of the Prandtl equation. Journal of the American Mathematical Society 23 (2), 591.
  * Gerard-Varet & Masmoudi (2013) Gerard-Varet, David & Masmoudi, Nader 2013 Well-posedness for the Prandtl system without analyticity or monotonicity. arXiv preprint arXiv:1305.0221 .
  * Goldstein (1948) Goldstein, S. 1948 On laminar boundary-layer flow near a position of separation. Q. J. Mech. Appl. Math. 1 (1), 43.
  * Grenier (2000) Grenier, Emmanuel 2000 On the nonlinear instability of Euler and Prandtl equations. Comm. Pure Appl. Math. 53, 1067–1091.
  * Grenier et al. (2014a) Grenier, Emmanuel, Guo, Yan & Nguyen, Toan T 2014a Spectral instability of characteristic boundary layer flows. arXiv preprint arXiv:1406.3862 .
  * Grenier et al. (2014b) Grenier, Emmanuel, Guo, Yan & Nguyen, Toan T 2014b Spectral stability of Prandtl boundary layers: an overview. arXiv preprint arXiv:1406.4452 .
  * Grenier et al. (2015) Grenier, E., Guo, Y. and Nguyen, T.T. 2015\. Spectral stability of Prandtl boundary layers: an overview. Analysis, 35(4), 343–355.
  * Gresho (1991) Gresho, P.M. 1991 Incompressible fluid dynamics: some fundamental formulation issues. Annual Review of Fluid Mechanics 23 (1), 413–453.
  * Heisenberg (1924) Heisenberg, W. 1924 Über Stabilität und Turbulenz von Flüssigkeitsströmen. Annalen der Physik 74, 577–627.
  * Hughes & Reid (1965) Hughes, TH & Reid, WH 1965 The stability of laminar boundary layers at separation. Journal of Fluid Mechanics 23 (04), 737–747.
  * Ingham (1984) Ingham, D.B. 1984 Unsteady separation. Journal of Computational Physics 53 (1), 90–99.
  * von Kármán (1921) von Kármán, T. 1921 Über laminare und turbulente Reibung. Z. ang. Math. Mech. 1 (4), 233–252, english translation in NACA TM 1092 (1946).
  * Kato (1972) Kato, T. 1972 Nonstationary flows of viscous and ideal fluids in $\mathbb{{R}}^{3}$. J. Funct. Anal. 9, 296–305.
  * Kato (1984) Kato, T. 1984 Remarks on zero viscosity limit for nonstationary Navier–Stokes flows with boundary. In Seminar on nonlinear partial differential equations, 85–98. MSRI, Berkeley.
  * Kelliher (2007) Kelliher, J.P. 2007 On Kato’s conditions for vanishing viscosity. Indiana University Mathematics Journal 56 (4), 1711–1722.
  * Kramer et al. (2007) Kramer, W., Clercx, H. & van Heijst, G. J. 2007 Vorticity dynamics of a dipole colliding with a no-slip wall. Phys. Fluids 19, 126603.
  * Lele (1992) Lele, S.K. 1992 Compact finite difference schemes with spectral-like resolution. Journal of Computational Physics 103 (1), 16–42.
  * Li et al. (1999) Li, X.S., Demmel, J.W., Gilbert, J.R., iL. Grigori, Shao, M. & Yamazaki, I. 1999 SuperLU Users’ Guide. Tech. Rep. LBNL-44289. Lawrence Berkeley National Laboratory, http://crd.lbl.gov/~xiaoye/SuperLU/. Last update: August 2011.
  * Lin (1967) Lin, C. C. 1967 The theory of hydrodynamic stability. Cambridge Univ Press.
  * Maekawa (2013) Maekawa, Y. 2013 Solution formula for the vorticity equations in the half-plane with application to high vorticity creation at zero viscosity limit. Advances in Differential Equations 18 (1/2), 101–146.
  * Maekawa (2014) Maekawa, Y. 2014 On the inviscid limit problem of the vorticity equations for viscous incompressible flows in the half-plane. Communications on Pure and Applied Mathematics 67(7), 1045–1128.
  * Nguyen van yen et al. (2009) 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, 89–107.
  * Nguyen van yen et al. (2011) 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.
  * Nguyen van yen et al. (2014) 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. Num. Math. 128, 301–338.
  * Obabko & Cassel (2002) Obabko, AV & Cassel, KW 2002 Navier–Stokes solutions of unsteady separation induced by a vortex. Journal of Fluid Mechanics 465, 99–130.
  * Oleinik (1966) Oleinik, O.A. 1966 On the mathematical theory of boundary layer for an unsteady flow of incompressible fluid. Journal of Applied Mathematics and Mechanics 30 (5), 951 – 974.
  * Orlandi (1990) Orlandi, Paolo 1990 Vortex dipole rebound from a wall. Phys. Fluids A 2, 1429–1436.
  * Peridier et al. (1991a) Peridier, V.J., Smith, FT & Walker, JDA 1991a Vortex-induced boundary-layer separation. part 1. the unsteady limit problem $\mathrm{Re}\to\infty$. J. Fluid Mech 232, 99–131.
  * Peridier et al. (1991b) Peridier, V.J., Smith, FT & Walker, JD 1991b Vortex-induced boundary-layer separation. part 2. unsteady interacting boundary-layer theory. Journal of Fluid Mechanics 232, 133–165.
  * Prandtl (1904) Prandtl, L. 1904 Über Flüssigkeitsbewegung bei sehr kleiner Reibung. In Proc. 3rd Inter. Math. Congr. Heidelberg, 484–491.
  * Prandtl (1921) Prandtl, L. 1921 Bemerkungen über die Entstehung der Turbulenz. Zeitschrift für Angewandte Mathematik und Mechanik 1 (6), 431–436.
  * Rayleigh (1880) Rayleigh, Lord 1880 On the stability, or instability, of certain fluid motions. Proc. London Math. Soc. 11 (1), 57–70.
  * Reed et al. (1996) Reed, Helen L, Saric, William S & Arnal, Daniel 1996 Linear stability theory applied to boundary layers. Annual Review of Fluid Mechanics 28 (1), 389–428.
  * Robinson (1991) Robinson, Stephen K. 1991 Coherent motions in the turbulent boundary layer. Annu. Rev. Fluid Mech 23, 601–639.
  * le Rond d’Alembert (1768) le Rond d’Alembert, Jean 1768 Paradoxe proposé aux géomètres sur la résistance des fluides. In Opuscules Mathématiques, , vol. 5, chap. XXXIV, 132–138. Briasson.
  * Sammartino & Caflisch (1998) Sammartino, Marco & Caflisch, Russel E. 1998 Zero viscosity limit for analytic solutions of the Navier–Stokes equation on a half-space. Comm. Math. Phys. 192, 433–491.
  * Schlichting (1979) Schlichting, H. 1979 Boundary layer theory. McGraw-Hill, New York.
  * Schubauer & Skramstad (1947) Schubauer, Galen Brandt & Skramstad, Harold Kenneth 1947 Laminar boundary layer oscillations and transitions on a flat plate. J. Aero. Sci. 14, 69–76, also NACA Report 909 (1948).
  * Sears & Telionis (1975) Sears, W. R. & Telionis, D. P. 1975 Boundary-layer separation in unsteady flow. SIAM J. Appl. Math. 28 (1), 215–235.
  * Smith (1988) Smith, FT 1988 Finite-time break-up can occur in any unsteady interacting boundary layer. Mathematika 35 (02), 256–273.
  * Sommerfeld (1909) 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, p. 116–124.
  * Stewartson (1974) Stewartson, K 1974 Multistructured boundary layers on flat plates and related bodies. Adv. Appl. Mech, vol. 14, 145–239. Academic Press.
  * Sutherland et al. (2013) Sutherland, D., Macaskill, C. & Dritschel, D. G. 2013 The effect of slip length on vortex rebound from a rigid boundary. Physics of Fluids (1994-present) 25 (9), 093104.
  * Temam & Wang (1997) Temam, R. & Wang, X. 1997 The convergence of the solutions of the Navier-Stokes equations to that of the Euler equations. Applied Mathematics Letters 10 (5), 29–33.
  * Tietjens (1925) Tietjens, Oskar 1925 Beiträge zur Entstehung der Turbulenz. Z. ang. Math. Mech. 5 (3), 200–217.
  * Tollmien (1929) Tollmien, W. 1929 Über die Entstehung der Turbulenz. Nachrichten von der Gesellschaft der Wissenschaften zu Göttingen, Mathematisch-Physikalische Klasse, 21–44, engl. trans. in NACA Tech. Rep. 609 (1931).
  * Van Dommelen (1981) Van Dommelen, Leon 1981 Unsteady boundary layer separation. PhD thesis, Cornell University.
  * Van Dommelen & Cowley (1990) 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, 593–626.
  * Van Dommelen & Shen (1980) Van Dommelen, L.L. & Shen, S.F. 1980 The spontaneous generation of the singularity in a separating laminar boundary layer. J. Comp. Phys. 38 (2), 125–140.