1 Introduction
In the framework of linear stability analysis, time-periodic flows can be classified into oscillatory flows (with zero mean) and pulsatile flows (with non-zero mean). In the latter case, the stability analysis of frozen base flow leads to accurate results, provided that the instability grows on a time scale much smaller than the period of the base flow. Therefore, the modulation of the mean flow can be considered as a secondary effect since the oscillatory part is small compared with the mean. However, for strict oscillatory cases, the time variation of the base flow is a fundamental effect and the quasisteady approximation becomes irrelevant.
The Stokes boundary layer is one of the simplest oscillatory shear flows. A sinusoidal pressure gradient, or a sinusoidal wall motion, is imposed on a fluid at rest in a semi-infinite domain. The three relevant physical parameters are the kinematic viscosity
${\it\nu}$
, the period of oscillations
$T$
and the amplitude of velocity variations
$U_{0}$
. The control parameter, namely the Reynolds number, takes the form
$$\begin{eqnarray}\mathit{Re}=U_{0}\sqrt{T/{\it\nu}},\end{eqnarray}$$
which represents the ratio of transient inertial force and viscous shear force.
Von Kerczek & Davis (Reference Von Kerczek and Davis1974) and Hall (Reference Hall1978) performed the first modal analysis, based on Floquet theory, but due to limited computer capability it was not possible to reach high enough Reynolds number values to find an instability. Some decades later Blennerhassett & Bassom (Reference Blennerhassett and Bassom2002) found the critical Reynolds number,
$\mathit{Re}_{c}=2511$
. However, that critical value is well above the thresholds observed in experiments; these discrepancies are thoroughly reviewed and discussed in Blennerhassett & Bassom (Reference Blennerhassett and Bassom2008). Ozdemir, Hsu & Balachandar (Reference Ozdemir, Hsu and Balachandar2014) have performed direct numerical simulations, and they confirm the observation of turbulent flows for Reynolds numbers above
$\mathit{Re}_{T}=900$
. A subcritical transition to turbulence is typical for wall-bounded shear flows, since this property is related to the non-normal nature of the linearised Navier–Stokes equations (Farrell Reference Farrell1988; Butler & Farrell Reference Butler and Farrell1992; Trefethen et al.
Reference Trefethen, Trefethen, Reddy and Driscoll1993; Schmid & Henningson Reference Schmid and Henningson1999; Schmid Reference Schmid2007). Consequently, the idealised analytical flow has quite different stability properties from those seen in practice. For example, the perturbation growth can be amplified by wall roughness (Blondeaux & Vittori Reference Blondeaux and Vittori1994; Vittori & Verzicco Reference Vittori and Verzicco1998) or a base flow variation (Thomas et al.
Reference Thomas, Blennerhassett, Bassom and Davies2015).
While the Floquet analysis (Blennerhassett & Bassom Reference Blennerhassett and Bassom2002) focuses on the asymptotic behaviour determined by the least stable mode, a general small disturbance is in fact a weighted combination of all eigenmodes. Moreover, because of the non-normality of the linearised stability operator, there is the potential for very large transient amplification of the disturbance energy, even in nominally stable flow conditions (Akhavan, Kamm & Shapiro Reference Akhavan, Kamm and Shapiro1991). The kinematic mechanism for the initial amplification of disturbance travelling with the shear flow was initially discovered by Orr (Reference Orr1907). In order to understand how an asymptotically stable flow becomes turbulent, Orr considered the stability analysis as an initial value problem. The result is a possibility of linear transient growth, eventually followed by an asymptotic (modal) decay. It is inferred that an initial growth could lead to disturbances so large as to trigger a nonlinear self-sustaining process, thereby triggering the turbulence. Nonetheless, the short-time dynamics raises the difficulty of determining the initial condition. As a solution, Farrell (Reference Farrell1988) computed the optimal perturbation, namely the initial condition that maximises the energy gain.
The present paper aims to apply this strategy to the subcritical transition in Stokes boundary layer flows. The optimisation method was developed by Andersson, Berggren & Henningson (Reference Andersson, Berggren and Henningson1999), Corbett & Bottaro (Reference Corbett and Bottaro2000) and Luchini (Reference Luchini2000), and is briefly recalled in § 2 for consistency; the results are presented in § 3. A critical discussion and some possible extensions are proposed in the conclusion.
2 Method
The Navier–Stokes equations are made dimensionless with the following characteristic scales: the period of oscillations
$T$
, the amplitude of the velocity
$U_{0}$
and the diffusive length scale
${\it\delta}=\sqrt{{\it\nu}T}$
. The flow is described by Cartesian coordinates (
$x,y,z$
), respectively the streamwise, wall-normal and spanwise directions. The configuration allows for two analytical solutions (Schlichting Reference Schlichting1979).
In the case of an oscillating wall, the solution (obtained by Stokes) is
$$\begin{eqnarray}U(y,t)=\cos (2{\rm\pi}t-\sqrt{{\rm\pi}}y)\text{e}^{-\sqrt{{\rm\pi}}y}.\end{eqnarray}$$
In the case of an oscillating pressure gradient, the solution is
$$\begin{eqnarray}U(y,t)=\cos (2{\rm\pi}t)-\cos (2{\rm\pi}t-\sqrt{{\rm\pi}}y)\text{e}^{-\sqrt{{\rm\pi}}y}.\end{eqnarray}$$
These two solutions ((2.1) and (2.2)) are shown in figure 1.
Figure 1. Laminar solution during half a period, for an oscillating wall (a) and oscillating external flow (b).
The Navier–Stokes equations, linearised around the laminar state
$U(y,t)$
, are
$$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\boldsymbol{u}=0,\\ \boldsymbol{u}_{t}+\mathit{Re}(\boldsymbol{U}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}\boldsymbol{u}+\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}\boldsymbol{U})=-\boldsymbol{{\rm\nabla}}p+{\rm\nabla}^{2}\boldsymbol{u}.\end{array}\right\}\end{eqnarray}$$
The Reynolds number is defined by (1.1). Perturbations (
$\boldsymbol{u}=[u,v,w]$
) are decomposed into the sum of Fourier modes in homogeneous directions, and the linear framework allows mode-by-mode analysis,
$\boldsymbol{u}(y,t)\text{e}^{\text{i}(k_{x}x+k_{z}z)}$
,
$k_{x}$
and
$k_{z}$
being respectively the streamwise and spanwise wavenumbers. Equations (2.3) are associated with no-slip boundary conditions at the wall (
$y=0$
),
$$\begin{eqnarray}\boldsymbol{u}=0,\end{eqnarray}$$
and an irrotational flow condition in the far field (
$y\rightarrow \infty$
),
$$\begin{eqnarray}\frac{\partial \boldsymbol{u}}{\partial y}=-\sqrt{k_{x}^{2}+k_{z}^{2}}~\boldsymbol{u}.\end{eqnarray}$$
The non-dimensional kinetic energy is such that
$$\begin{eqnarray}E(t)=\int _{y}\frac{\boldsymbol{u}^{\ast }\boldsymbol{\cdot }\boldsymbol{u}}{2}\,\text{d}y,\end{eqnarray}$$
with
$\ast$
denoting complex conjugate.
The energy equation for the perturbations is obtained from the integrated momentum equation,
$$\begin{eqnarray}\frac{\text{d}E}{\text{d}t}=-\mathit{Re}\int _{y=0}^{\infty }\overline{uv}\frac{\partial U}{\partial y}\,\text{d}y-\int _{y=0}^{\infty }u_{i,j}u_{i,j}\,\text{d}y.\end{eqnarray}$$
The first right-hand-side term is the production term, where
$\overline{uv}$
represents the average along the homogeneous directions (
$x,z$
) of the product of the streamwise and wall-normal velocities, and
$\partial U/\partial y$
denotes the laminar shear. The last term, always negative, designates the dissipation of kinetic energy by the viscous effect.
From (2.7), it can be stated that the two Stokes solutions ((2.1) and (2.2)) lead to identical energy growth. In particular, this is true for the less stable mode, as shown, with a different argument, by Blennerhassett & Bassom (Reference Blennerhassett and Bassom2002): they proved that eigenvalues are invariant under the transformation
$U(y,t)\rightarrow U(y,t)-\cos (2{\rm\pi}t)$
.
Given that the basic state can grow (or decay) simultaneously with the growth of a disturbance, it is sometimes more appropriate to consider a relative criterion for instability defined as the energy ratio of the perturbation to the base flow (Conrad & Criminale Reference Conrad and Criminale1965). This is especially true for the case where the kinetic energy of the base flow would become null, so that the infinitesimal assumption for the perturbation would become questionable. However, this consideration is not significant in the present case because the energy of the laminar state is strictly positive.
The concept of optimal perturbations was first introduced by Farrell (Reference Farrell1988) and further developed by Butler & Farrell (Reference Butler and Farrell1992). The method was initially based on the computation of the full spectrum of the linearised Navier–Stokes operator. A significant improvement was proposed by Andersson et al. (Reference Andersson, Berggren and Henningson1999), Corbett & Bottaro (Reference Corbett and Bottaro2000) and Luchini (Reference Luchini2000), with an optimisation based on the continuous adjoint operator. This method, initially used for steady base flow, can easily be applied to base flows with arbitrary time dependence (Zhao, Ghidaoui & Kolyshkin Reference Zhao, Ghidaoui and Kolyshkin2007). The applications of adjoint methods for linear stability analysis are reviewed in depth by Luchini & Bottaro (Reference Luchini and Bottaro2014).
If we formally denote by
$L\boldsymbol{u}=0$
the linearised equations (2.3), then the optimisation problem is to maximise the energy
$E(t)$
subject to the constraint
$L\boldsymbol{u}=0$
. The Lagrangian functional is based on an energy-like norm and reads
$$\begin{eqnarray}\mathscr{L}(\boldsymbol{u},\boldsymbol{u}^{\dagger })=E(t)+\int _{V}\boldsymbol{u}^{\dagger \ast }\boldsymbol{\cdot }(L\boldsymbol{u})\,\text{d}V,\end{eqnarray}$$
with
$\boldsymbol{u}^{\dagger }$
the Lagrange multiplier. The optimal state is reached for the stationary point for the Lagrange function, which gives us the adjoint problem:
$$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\boldsymbol{u}^{\dagger }=0,\\ -\boldsymbol{u}_{t}^{\dagger }-\mathit{Re}(\boldsymbol{U}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}\boldsymbol{u}^{\dagger }-\boldsymbol{u}^{\dagger }\boldsymbol{\cdot }(\boldsymbol{{\rm\nabla}}\boldsymbol{U})^{\text{T}})=-\boldsymbol{{\rm\nabla}}p^{\dagger }+{\rm\nabla}^{2}\boldsymbol{u}^{\dagger }.\end{array}\right\}\end{eqnarray}$$
The adjoint equations (2.9) are associated with the following boundary conditions: no-slip conditions at the wall,
$$\begin{eqnarray}\boldsymbol{u}^{\dagger }=0,\end{eqnarray}$$
and an irrotational flow condition in the far field,
$$\begin{eqnarray}\frac{\partial \boldsymbol{u}^{\dagger }}{\partial y}=-\sqrt{k_{x}^{2}+k_{z}^{2}}~\boldsymbol{u}^{\dagger }.\end{eqnarray}$$
Practically, to identify the flow state at
$t=T_{0}$
producing the largest disturbance growth at any given
$T_{f}>T_{0}$
, a repeated numerical integration of the direct and adjoint equations is used, coupled with transfer conditions. The direct problem (the linearised Navier–Stokes equations) is integrated from a given initial condition at
$t=T_{0}$
up to the final target time
$t=T_{f}$
; afterwards the adjoint problem is integrated backwards in time. The initial condition for the adjoint problem, at
$t=T_{f}$
, is
$\boldsymbol{u}^{\dagger }=\boldsymbol{u}$
, whereas the direct problem is initiated, at
$t=T_{0}$
, with
$\boldsymbol{u}=\boldsymbol{u}^{\dagger }$
. The perturbations are normalised such that the energy is equal to one at
$t=T_{0}$
. The direct–adjoint iterative procedure is stopped when the energy
$E(T_{f})$
has converged to a value, with a precision below
$10^{-8}$
.
The equations are solved with a Chebyshev spectral collocation method in the inhomogeneous direction (
$y$
). The Gauss–Lobatto–Chebyshev grid
${\it\xi}$
is expanded from the interval
$[-1,1]$
to the physical domain
$[0,y_{\infty }]$
through a linear mapping,
$y=y_{\infty }(1-{\it\xi})/2$
. The divergence-free condition is achieved with a prediction–projection scheme. The pressure is approximated by polynomials of two units lower order than for the velocity, in order to compute a pressure unpolluted by spurious modes with only one collocation grid (Botella Reference Botella1997). The time marching combines a fourth-order Adams–Bashforth scheme and a fourth-order backward differentiation scheme, with the viscous term treated implicitly (Ascher, Ruuth & Wetton Reference Ascher, Ruuth and Wetton1995). The numerical parameters are
$y_{\infty }=10$
,
$N_{y}=100$
and
${\rm\Delta}t=10^{-5}$
.
3 Results
As preliminary results, the neutral stability results are obtained with
$T_{0}=0$
,
$k_{z}=0$
and
$T_{f}=4$
. The critical Reynolds number is obtained when the final energy is close to the initial energy
$E(T_{f})=E(T_{0})$
, provided that
$T_{f}$
is sufficiently large. The computed values are
$Re_{c}=2511$
and
$k_{x}=0.667$
, in line with the values obtained by Blennerhassett & Bassom (Reference Blennerhassett and Bassom2002) and Thomas et al. (Reference Thomas, Bassom, Blennerhassett and Davies2010).
In the subcritical regime, the parametric study is realised by varying the five parameters: the streamwise (
$k_{x}$
) and spanwise (
$k_{z}$
) wavenumbers, the initial (
$T_{0}$
) and final (
$T_{f}$
) times, and the Reynolds number. For all cases, the optimal perturbation is two-dimensional; in other words, the spanwise velocity component and the spanwise wavenumber are both zero.
The optimal streamwise wavenumber is depicted in figure 2(a); the optimal value remains close to the critical value (
$k_{x}=0.667$
). For low Reynolds number values, the optimal time to inject the perturbation, namely
$T_{0}$
, does not exactly correspond to the beginning of the deceleration phase (
$t=0$
), see figure 2(b). Indeed, when the deceleration starts, the mean shear is primarily located close to the wall, where viscous damping is more effective; as a result, amplification is hampered. This results in a compromise for the optimal start of the Orr mechanism. For relatively large Reynolds numbers, the optimal part of the cycle, shown in figure 2(b), spreads throughout half a period: starting from the beginning of a deceleration phase up to the end of the next acceleration phase.
Figure 2. Optimal parameters for various Reynolds numbers: optimal streamwise wavenumber (a) and temporal interval (b). The optimal spanwise wavenumber is zero.
The maximum energy amplification increases exponentially with the Reynolds number (see figure 3). This scaling is characteristic of velocity profiles with reverse flow, for example the flow over a backward-facing step (Blackburn, Barkley & Sherwin Reference Blackburn, Barkley and Sherwin2008) and the flow through an expanding pipe (Cantwell, Barkley & Blackburn Reference Cantwell, Barkley and Blackburn2010). The transient growth of energy is fairly large, so that the modal path could be bypassed in supercritical regimes and replaced by this linear Orr mechanism, even for low disturbance levels.
Figure 3. Maximum transient energy growth for various Reynolds numbers. The perturbations are normalised such that
$E(t=T_{0})=1$
.
The present results differ significantly from those obtained for steady shear flows (Butler & Farrell Reference Butler and Farrell1992; Schmid & Henningson Reference Schmid and Henningson1999). In the latter case, streamwise vortices (with
$k_{x}=0$
and
$k_{z}\neq 0$
) transform into streaks. The maximum disturbance energy, mostly carried by the streaks, scales as the square of the Reynolds number increases, which is much lower than the scale obtained for oscillatory flows. Moreover, for steady flows, the effect of acceleration is to damp growth (Corbett & Bottaro Reference Corbett and Bottaro2000); in the present case the perturbations seem to be less sensitive to this damping effect.
The underlying mechanism for the transient amplification in oscillatory flows corresponds to an Orr mechanism, as can be seen in figure 4.
Figure 4. Contours of the spanwise vorticity of the perturbations and the corresponding base flow profile for
$t=T_{0}$
(a,b) and
$t=T_{f}$
(c,d), at
$Re=1000$
.
The mean shear is initially tilted against the direction of the initial disturbances. In this configuration, the structure borrows energy from the mean flow via the Reynolds stress production term (see (2.7)). The inclination of the spanwise vortices (figure 4) illustrates that the term
$\overline{uv}$
is negative as long as the shear
$\partial U/\partial y$
is positive. After half a period, the mean flow is tilted and aligned with the vorticity of the perturbation; as a consequence, the energy of the disturbances decreases. The transient growth mechanism is presented for the oscillating pressure gradient case, but the same mechanism is observed for the case of an oscillating wall.
In order to complete the qualitative results, values are presented in table 1, in particular for the
$\mathit{Re}=1000$
case, which corresponds to the observed value for transition to turbulence.
Table 1. Optimal parameters for three Reynolds number values.
The convergence of the present numerical method has been assessed by comparing the maximum energy growth computed with refined numerical parameters:
$y_{\infty }=15$
,
$N_{y}=200$
and
${\rm\Delta}t=10^{-6}$
. For example, at
$\mathit{Re}=1000$
, the obtained value is
$E(T_{f})=3.7741\times 10^{6}$
. Comparison with the value in table 1 gives an estimate for the convergence to four significant digits.
Nonlinear simulations have been conducted at
$\mathit{Re}=1000$
and 1200, with the optimal parameters given in table 1. In the following, the flow is driven by an oscillating pressure gradient, but the results are qualitatively similar to those obtained with an oscillating wall (not shown). The streamwise and spanwise lengths of the domain are
$2{\rm\pi}/k_{x}$
, with
$k_{x}$
the optimal streamwise wavenumber. The numerical parameters are
$y_{\infty }=10$
,
$N_{y}=100$
,
${\rm\Delta}t=10^{-5}$
and
$N_{x}=N_{z}=128$
Fourier modes in the homogeneous directions. These numerical parameters are consistent with those used by Ozdemir et al. (Reference Ozdemir, Hsu and Balachandar2014). The initial states consist of the unsteady base flow profile plus the optimal perturbation, normalised with the prescribed energy
$E_{0}=10^{-7}$
and
$E_{0}=10^{-9}$
for
$\mathit{Re}=1000$
and 1200 respectively. A white Gaussian divergence-free noise is superimposed in order to break the spanwise invariance, with a prescribed energy
$E_{1}=10^{-3}E_{0}$
. The two simulations similarly lead to turbulent flows; the case
$\mathit{Re}=1000$
is illustrated in figure 5.
Figure 5. The temporal evolution of the kinetic energy of the perturbations, at
$\mathit{Re}=1000$
, computed in three different frameworks, namely the three-dimensional nonlinear simulation (DNS 3D), the two-dimensional nonlinear simulation (DNS 2D) and the linear evolution.
In figure 5, the evolution energy of the perturbations obtained from the three-dimensional simulation is compared with the two-dimensional case and the linear evolution. The nonlinear evolution follows the linear curve over four decades, at which point the base flow is slightly distorted by quadratic interaction. As a consequence, the amplification saturates and the maximal energy gain does not reach the linear value. While the two-dimensional simulation shows a relaminarisation, with a slow decay of the mean flow defect, the three-dimensional case shows that the two-dimensional wave destabilises into three-dimensional disturbances, with an evolution towards turbulence. As for steady flows, three-dimensional fluctuations are necessary to observe a self-sustaining wall turbulence in subcritical regimes. Indeed, as a common feature of the near-wall turbulence, the development of streamwise low-speed streaks is observed at the end of the transition to turbulence (Jimènez & Moin Reference Jimènez and Moin1991; Hamilton, Kim & Waleffe Reference Hamilton, Kim and Waleffe1995). The amplification is fast; the whole transition can be realised within one deceleration phase. Following the transient, the flow is close to the laminar state in the accelerated phase, but turbulent bursts are triggered during the decelerated phase of the cycle. This transition towards turbulence is similar to the simulations initialised with non-optimal perturbations by Ozdemir et al. (Reference Ozdemir, Hsu and Balachandar2014).
4 Conclusion
The subcritical instability of oscillatory Stokes flows was investigated with a linear optimal perturbation analysis. The parametric study indicates that the disturbances are able to efficiently extract energy from the mean shear by transporting momentum down the mean velocity gradient. The results show a strong Orr mechanism, with amplifications increasing exponentially with the Reynolds number. The optimal disturbances take the form of two-dimensional spanwise vortices, with streamwise wavenumber close to the value
$k_{x}=0.7$
.
The optimisation strategy completes the modal (asymptotic) analysis by Blennerhassett & Bassom (Reference Blennerhassett and Bassom2002), in particular for subcritical Reynolds numbers. As the Orr mechanism is non-modal by nature, this process can only be observed by considering an initial value problem. Nonetheless, an optimisation strategy for computing the initial condition raises the issue of the frequency of occurrence of such an event. While the unstable mode eventually emerges in supercritical conditions, the transient amplification results in a subtle combination of eigenmodes, determined by the receptivity process.
The present results, applied to classical oscillatory boundary layer flows, can be extended to more complex flows by using a linear combination of the Stokes solutions. For example, Spalart (Reference Spalart1989) used Stokes solutions for the velocity components in the streamwise and spanwise directions in order to create simple turbulent boundary layer flows with three-dimensional statistics. A possible extension of the research presented in the present paper would be the stability analysis of the superposition of two Stokes solutions, with an irrational frequency ratio.
