Revisiting the momentary stability analysis of the Stokes boundary layer

Abstract The stability of the boundary layer generated by the harmonic oscillations of a plate in its own plane in a fluid otherwise at rest (Stokes boundary layer) is investigated by considering the time development of perturbations of small amplitude and introducing a momentary criterion of instability. The temporal scale of the perturbations is assumed to be much smaller than the period of the plate oscillations because transition takes place at values of the Reynolds number much larger than one. The results confirm that the Stokes boundary layer is linearly unstable when the Reynolds number is larger than a first critical value $R_{\delta ,c1}$ equal to $85$, which is almost coincident with that determined by Von Kerczek & Davis (J. Fluid Mech., vol. 62, issue 4, 1974, 753–773) for a Stokes boundary layer in a fluid domain which is bounded by a second fixed plate at a distance $h^{*}$ from the oscillating one. For values of the Reynolds number close to $R_{\delta ,c1}$, the instability is restricted to phases close to the inversion of the plate velocity. When the Reynolds number becomes larger than a second threshold value $R_{\delta ,c2}$ close to $200$, the instability rapidly pervades a large part of the cycle. However, only when the Reynolds number becomes larger than a third critical value $R_{\delta ,c3}$ equal to $850$, is the instability present during the whole cycle. Heuristically, these three critical values of the Reynolds number can be associated with the transition from the laminar regime to the disturbed laminar, the intermittently turbulent and fully turbulent regimes.

the harmonic oscillations of a plate in an otherwise fluid at rest (Stokes boundary layer) were carried out and continue to be carried out. Indeed, the Stokes boundary layer is the prototype of the unsteady boundary layers which are present in a large number of applications ranging from coastal engineering to biomedical sciences.
The flow in the laminar regime is well known and it was determined by Stokes (1851). More recent studies address the problem of detecting transition to turbulence and studying turbulence characteristics.
Two different approaches are commonly used to investigate the stability of the Stokes boundary layer and transition to turbulence.
The first approach is based on Floquet theory and it looks for the net growth/decay of a perturbation of the basic flow over the entire period of oscillation.
The second approach, also known as the quasi-steady approach, assumes that the perturbation growth/decay takes place on a time scale much smaller than the period of the basic flow. Hence, the eigenvalue problem describing the time development of the perturbation is solved at each phase of the oscillation cycle to obtain an instantaneous growth rate. This approach, which is justified only when the Reynolds number of the basic flow is assumed to have asymptotically large values, introduces a momentary criterion of instability and defines the flow to be unstable when a phase of the cycle exists at which the perturbation grows (Cowley 1987;Hall 2003).
The approach which looks at the phase averaged growth of the perturbations was first used by Hall (1978). Hall (1978) solved the problem following the approach of Seminara & Hall (1976) but, because of the limited power of the computers at that time, he could investigate only moderate values of the Reynolds number and he found no instability of the Stokes boundary layer up to values of the Reynolds number equal to 800. Hereinafter, the Reynolds number R δ is defined by where U * 0 is the amplitude of the velocity oscillations of the plate and δ * = √ 2ν * /ω * is a conventional viscous length defined using the kinematic viscosity ν * of the fluid and the angular frequency ω * of the plate oscillations. Later, using the same approach employed by Hall (1978), Blennerhassett & Bassom (2002) could determine the behaviour of the perturbations for larger values of the Reynolds number and they found unstable modes when R δ is larger than 1416. Incidentally, let us point out that Blennerhassett & Bassom (2002) defined the Reynolds number as U * 0 / √ 2ν * ω * , a value which is half of R δ , and other authors use the Reynolds number Re = U * 2 0 /ν * ω * = R 2 δ /2. The approach based on the momentary criterion of stability was first used by Blondeaux & Seminara (1979), who predicted the growth of harmonic components characterized by a wavelength equal to approximately 12.5δ * when the Reynolds number becomes larger than 85 (since the analysis developed by Blondeaux & Seminara (1979) and the results they obtained are described in Italian, a brief summary of the analysis of Blondeaux & Seminara (1979) is provided in the following).
At this stage, it is worth pointing out that the analyses of Hall (1978), Blennerhassett & Bassom (2002) and Blondeaux & Seminara (1979) consider the 'infinite' case, i.e. the boundary layer generated by a flat plate oscillating in a semi-infinite domain. As discussed in Hall (1978), the results provided by the stability analyses of the finite case (the oscillatory flow in a two-dimensional duct) differ from those of the infinite case.
These linear stability analyses attempt to explain the experimental observations that show the existence of different 'flow regimes'. Even though in the scientific literature different terms are used and different definitions are given, we think it is reasonable to identify (i) the laminar regime, (ii) the disturbed laminar regime, (iii) the intermittently turbulent regime, (iv) the fully turbulent regime.
The disturbed laminar regime occurs when the Reynolds number exceeds a first critical value R δ,c1 (R δ,c1 100) but it stays smaller than a second critical value R δ,c2 and it is characterized by the transient appearance of small amplitude perturbations of the laminar flow. However, the overall flow does not differ from that of the laminar regime (see the 'weakly turbulent regime' defined by Hino, Sawamoto & Takasu (1976)).
The intermittently turbulent regime takes place when the Reynolds number is larger than R δ,c2 and is characterized by the presence of violent bursts of turbulence which appear during the decelerating phases of the cycle. Experimental observations show that R δ,c2 falls between 500 and 600. Then, the flow recovers a laminar-like behaviour during the accelerating phases (see the 'conditional turbulent regime' defined by Hino et al. (1976)).
The duration of the turbulence bursts increases as the Reynolds number is increased. In particular, the measurements of Hino et al. (1983), Jensen, Sumer & Fredsøe (1989 and Akhavan, Kamm & Shapiro (1991a) show that turbulence is present throughout the whole cycle when the Reynolds number R δ is larger than third a critical value R δ,c3 which falls close to 750. The reader can look at the measurement of the turbulent kinetic energy plotted in figures 11-13 of Hino et al. (1983), figures 19-20 of Jensen et al. (1989) and figures 9-13, 24-26 of Akhavan et al. (1991a).
In a few papers it is stated that so far no experimental observation exists that shows turbulence presence throughout the whole oscillatory cycle. However, this statement is due to a misinterpretation of one of the conclusions written by Jensen et al. (1989) which reads '. . . the present experiments indicate that, even at Re as large as 1.6 × 10 6 , there is still some portion of the half-cycle (namely ωt ≤ 45 • ) where the flow regime is not a fully developed turbulent one'. With the term 'fully developed turbulent regime', Jensen et al. (1989) means a flow regime characterized by a turbulence level so high that the velocity profile close to the wall follows the log-law; they do not mean the simple presence of turbulence that was significant during the whole cycle. Indeed, figures 19 and 20 of Jensen et al. (1989), where the results of Hino et al. (1983) are also plotted, show that the root mean square values of the velocity oscillations assume significant values during the whole cycle for Re equal to 2.8 × 10 5 and 5.0 × 10 5 , i.e. R δ = 748 and R δ = 1000, respectively. Turbulence presence throughout the oscillatory cycle was confirmed by one of the authors of the paper of Jensen et al. (1989) (B.M. Sumer (2020, private communication).
A quantitative comparison of the experimental observations with the results of the stability analyses available in the scientific literature is not entirely satisfactory. Even though the momentary criterion of stability can explain the appearance of perturbations when the Reynolds number becomes larger than the first critical value R δ,c1 , and the theoretical value R δ,c1 = 85 is not too far from the experimental value which ranges around 100, the theoretical analyses do not explain the cyclic appearance of the turbulence bursts for values of R δ larger than a value falling between 500 and 600. In fact, according to Floquet theory, any initial perturbation should decay after a few cycles up to values of R δ equal to 1416 (Blennerhassett & Bassom 2002). A possible mechanism able to trigger the periodic appearance of turbulence was pointed out by Blondeaux & Vittori (1994) who showed that turbulence appearance can be triggered by wall imperfections through a resonance phenomenon. Blondeaux & Vittori (1994) found that when the Reynolds number is larger than a threshold value ranging around 100, an instant within the decelerating phases of the cycle does exist such that a waviness of infinitesimal amplitude can induce large amplitude perturbations of the Stokes flow if its wavelength has particular values. Later, further investigations of this resonant mechanism was made by Luo & Wu (2010) and Scandura (2013).
In this paper, we revisit the stability analysis of the Stokes flow made by Blondeaux & Seminara (1979) using a momentary criterion of instability and we evaluate what are the portions of the cycle during which the Stokes flow turns out to be unstable. The results of the analysis allow us to determine not only the critical value of the Reynolds number above which the Stokes boundary layer turns out to be momentarily unstable but also other critical values above which the intermittently turbulent and fully turbulent regimes are expected to be observed.

Formulation of the problem and solution approach
Let us consider a flat plate which bounds a viscous fluid and oscillates in its own plane with velocity U * 0 cos(ω * t * ), where U * 0 and ω * denote the amplitude and the angular frequency of the velocity oscillations, and t * is time. Moreover, let us introduce the viscous length δ * = √ 2ν * /ω, where ν * is the kinematic viscosity of the fluid. Let us introduce a coordinate system (x * , y * , z * ) with the x * -and y * -axes lying on the plate, the x * -axis being in the direction of the plate oscillations, and the z * -axis being vertical and pointing upwards.
The hydrodynamic problem is posed by continuity and Navier-Stokes equations along with the no-slip condition at the plate surface and the boundary condition that enforces the vanishing of the velocity far from the plate. Let us introduce dimensionless variables such that the spatial coordinates are scaled with the viscous length δ * and time is scaled using the angular frequency of the fluid oscillations, Moreover, let us denote with (u, v, w) and p the dimensionless velocity components and the pressure field, respectively, defined by where ρ * is the density of the fluid. Then, the hydrodynamic problem is posed by where the force per unit mass is assumed to be conservative and the modified pressure p * m is introduced such that p * m = p * − ρ * ϕ * , ϕ * = −g * z * being the dimensional potential function of the body force and g * being the gravitational acceleration.
It can be easily verified that flow field induced by the plate oscillations is described by where p m,0 is the leading-order term of the modified pressure, i.e. the actual pressure minus the hydrostatic contribution.
To investigate the stability of this basic flow, let us consider a perturbation of small amplitude superimposed on the basic flow and let us determine its time development. Since, according to Squire's theorem, two-dimensional perturbations of the flow field are more unstable than three-dimensional perturbations, let us assume (u, v, w, p m where is assumed to be much smaller than unity. Neglecting terms of order 2 , continuity and momentum equations read The continuity equation is satisfied if a stream function ψ is introduced such that (2.12) Then, eliminating the pressure field from (2.10) and (2.11), the linearized vorticity equation is obtained as (2.13) Because of the linearity of (2.13), it is possible to consider the generic streamwise harmonic components of the perturbation. Moreover, because the growth of perturbations of the basic laminar flow is observed when the Reynolds number is much larger than unity, the amplitude of the flow perturbations can be assumed to grow on a temporal scale τ much faster than that characterizing the time development of the basic flow. Hence, a 'momentary' criterion for the instability of the basic flow can be used (Shen 1961). Indeed, the Reynolds number can be thought to be the ratio between the characteristic temporal scale 1/ω * of the fluid oscillations and the convective temporal scale δ * /U * 0 of the time development of the perturbation. Hence, let us consider a solution of (2.13) in the form where h.o.t. indicates higher-order terms. In (2.14), α is the wavenumber of the generic Fourier component of the flow perturbation in the streamwise direction, the real part c r of the complex quantity c is its wave speed and the imaginary part c i controls its growth/decay rate. At the leading order of approximation, when (2.14) is plugged into (2.13), the following homogeneous differential equation is obtained: where the operator N 2 is defined by N 2 = ∂ 2 /∂z 2 − α 2 . The boundary conditions force the vanishing of the velocity at the wall and far from it At this stage, it is worth pointing out that the viscous term is retained in (2.15) even if it is of order 1/R δ , because it turns out to be significant, at the leading order of approximation, within a viscous layer close to the wall and within critical layers. As pointed out by Blondeaux & Vittori (1994), a formal approach would require the solution of the inviscid version of (2.15) and, then, the matching of the inviscid solution with the solution which holds in the viscous and critical layers. However, such a procedure would involve much tedious and heavy algebra which can be avoided by the direct solution of (2.15). Blondeaux & Vittori (1994) showed that the direct solution of (2.15) coincides, to the required order of approximation, with the solution obtained on the basis of the rigorous matched asymptotic approach, previously briefly outlined.
A non-vanishing solution of the hydrodynamic homogeneous problem posed by (2.15) and its boundary conditions can be found by forcing an eigenrelation which provides the value of c as a function of t, α and R δ .
To find the eigenvalues c and the corresponding eigenfunctions f , the approach employed by Blondeaux & Seminara (1979) is employed. The stream function ψ is written in the form where σ = α 2 − 2 iαR δ c and the coefficients a n,m and b n,m can be written as functions of a 0,0 and b 0,0 by means of the recursive relationships a n,m , b n,m = iαR δ μ 2 e it − (α, σ ) − C n,m 2 − α 2 − 2i a n−1,m , b n−1,m Moreover, μ = 0 when m = n, whereas μ = 1 when 0 ≤ m ≤ n − 1. Similarly, ζ = 0 when m = 0, and ζ = 1 when 1 ≤ m ≤ n.
The boundary conditions at the wall lead to n m=0 a n,m a 0,0 , A non-trivial solution of the system (2.20) can be found if and only if which provides the eigenrelation c = c(t, α, R δ ). This flow field vanishes when z tends to infinity unless c r = 0 and c i < −(αR δ /2). As discussed in Blondeaux & Seminara (1979), when c r = 0 and c i < −(αR δ /2), a continuous spectrum of eigenvalues exists but the eigenfunctions do no decay moving far away from the oscillating wall. Indeed, for such values of c the solution of the eigenvalue problem becomes f = a e −αz + b cos α 2 + αR δ c i z + c sin α 2 + αR δ c i z + exponential decaying terms.
(2.23) Therefore, the continuous spectrum will be no longer considered in the following.

The results
The eigenfunction f and the eigenvalues c are evaluated by truncating the summation (2.17) after N terms, with N large enough that the addition of further terms does not modify the results. The value of N depends on R δ and α. For example, for α = 0.5 and R δ = 100, already N = 15 provides reliable results (five significant digits), but for α = 0.5 and R δ = 1000, values of N larger than 35 are required to obtain the same accuracy. However, the evaluation of A 11 , A 12 , A 21 and A 22 requires a relatively small effort even when the accurate evaluation of ψ asks for a large number of terms in (2.17). It is worth pointing out that increasing values of N demand the evaluation of a n,m and b n,m with more digits. For example, for α = 0.5, R δ = 100 and t = 0.0, the single-precision floating-point format is enough, but for α = 0.5, R δ = 1500 and t = 0.0, the quadruple-precision is required.
Because of the small computational costs required by the numerical procedure used to (2) obtain the eigenvalues c and the eigenfunctions f , the results described in the following were obtained by fixing N = 75 and using the quadruple-precision for all the values of α and R δ . Before discussing in detail the quantitative results provided by the solution of the eigenvalue problem described in the previous section, let us point out that for any eigenvalue c(t) = c r (t) + ic i (t), there are other eigenvalues equal to minus the complex conjugate of c at phases of the cycle equal to t ± π (c r (t) + ic i (t) = −c r (t ± π) + ic i (t ± π)). Figure 1 shows the real and imaginary parts of c as a function of the phase t of the cycle for α = 0.5 and two different values of the Reynolds number, namely R δ = 80 and R δ = 90, which are just below and above the critical value determined by Blondeaux & Seminara (1979), namely R δ,c1 = 85. For the smaller value of R δ , the growth rate of the perturbation, i.e. the imaginary part of c, is always negative and the perturbation always decays. On the other hand, for the larger value of R δ , small parts of the cycle exist such that the perturbation grows.  The size of the unstable parts of the cycle increases as the Reynolds number is increased, as it appears looking at figure 2 where the eigenvalue c i is plotted versus t for α = 0.5 and R δ = 150 and R δ = 300. Figure 3 shows the eigenfunction f (see (2.17)) for R δ = 300, α = 0.5, t = 3π/2 and t = 1.8396. At t = 3π/2, the eigenfunction is significant in a layer, the thickness of which is similar to that of the Stokes boundary layer, i.e. the |f | is significant up to values z of order 1. At t = 1.8396, the discrete spectrum of eigenvalues is going to merge into the continuous spectrum and the eigenfunction tends to pervade the whole fluid domain since for smaller values of t no mode exists that decays moving away from the oscillating plate. By analysing at the velocity profiles of the basic flow at different phases of the cycle, it appears that the instability of the Stokes flow is predicted when an inflection point is present in a region of the flow where both the velocity and its derivative are large and the Reynolds number is large enough to lead to a significant interaction between the perturbations and the basic flow. An investigation of the parameter space showed that the Stokes boundary layer is 'momentarily' unstable during the oscillatory cycle for values of R δ larger than 85 and the most unstable wavenumber α c,1 is equal to 0.5. Figure 4 shows the maximum value of the imaginary part c i of c as a function of α and R δ irrespective of the phase of the cycle at which the maximum takes place. The line characterized by c i = 0 is the marginal stability curve in the plane (α, R δ ). Figure 5 shows the phase t ins of the cycle at which the instability first occurs as a function of the Reynolds number. Close to the critical conditions, the instability is predicted close to t = π/2 + nπ with n = 0, 1, 2, . . . but the instability takes place earlier as the Reynolds number is increased.
These results were already obtained by Blondeaux & Seminara (1979). However, Blondeaux & Seminara (1979) did not look at the behaviour of c for larger values of R δ . Figure 6 shows the imaginary part of c as a function of t for α = 0.5 and R δ = 500 and 1500. For the smallest value of R δ , there is a time interval during which even the most unstable mode is characterized by a negative value of c i , implying a momentary decay of the perturbation. On the other hand, for R δ = 1500, there is an unstable mode during the whole cycle. In figure 6, the points D indicate the phase of the cycle such that the perturbation which propagates in the positive direction of the x-axis has the same growth rate as the perturbation which propagates in the negative direction. To determine the critical value R δ,c3 which leads to the growth of perturbations at any phases of the cycle, the evaluation of c i is carried out for many values of R δ and α. Figure 7 shows the value c i at the point D (see panel (b) of figure 6) plotted in the plane (α, R δ ). Looking at figure 6, it can be realized that the instability pervades the whole oscillation cycle when c i at the point D is zero or positive. It follows that the results plotted in figure 7 allow the determination of the value R δ,c3 of the Reynolds number such that, for R δ larger than R δ,c3 , a growing perturbation is present during the whole cycle. It turns out that R δ,c3 = 818 and the most unstable perturbation is characterized by the wavenumber α c,3 = 0.4. Hence, it can be inferred that the fully turbulent regime is predicted by the present analysis when the Reynolds number is larger than 818.
In figure 8, the value of c i is plotted for α = 0.5 and different values of R δ . When the Reynolds number is not too far from R δ,c1 , the instability is restricted to small parts of the oscillatory cycle. However, when the Reynolds number is increased beyond 200, the unstable time interval rapidly increases, and for R δ 550 the instability pervades a large part of the decelerating phases and the whole accelerating phases. Indeed, the lower limit t 1 of the unstable time intervals does not change significantly when the Reynolds number is increased, whereas the upper limit t 2 moves quickly to 2π. We remind the reader that c i (t − π) = c i (t) is also an eigenvalue which is not plotted in figure 8 for the sake of clarity. Hence, the instability phases also range between t 1 ± nπ and t 2 ± nπ. The existence of a further critical value of the Reynolds number R δ,c2 can be inferred, beyond which, not only a growing perturbation exists, but transition to turbulence is promoted during a large time interval. The different qualitative behaviour of the perturbations, when the Reynolds number moves from 200 to 550 and beyond, clearly appears when the values of t 1 and t 2 , i.e. the boundaries of the unstable time intervals, are plotted versus the Reynolds number. Figure 9 clearly shows that, while t 1 decreases continuously with R δ , t 2 is characterized by a sudden and rapid growth when R δ becomes larger than approximately 200. The rapid growth of the unstable phases is more evident if t 2 and dt 2 /dR δ are plotted versus the Reynolds number. Indeed figure 10 shows that the value of dt 2 /dR δ decreases from large values to small values when R δ moves from 85 to 187.5 but it suddenly jumps to larger values when R δ moves from 187.5 to 200, and decays again for values of R δ larger than 200.

Conclusions
The present momentary stability analysis extends to the 'infinite' case the analysis of Von Kerczek & Davis (1974) who considered the 'finite' case, introducing a wall at a distance h from the oscillating wall which generates the Stokes boundary layer.
Even though the present results come from a linear analysis of the time development of perturbations of the Stokes boundary layer and, certainly, it cannot be claimed that they describe turbulence dynamics, the present findings explain the presence of the different flow regimes observed in the laboratory experiments (see, amongst others, Hino et al. (1976), Hino et al. (1983), Jensen et al. (1989) and Akhavan et al. (1991a)).
Indeed, the present results show the following.
(i) Perturbations of the laminar flow start to appear when R δ becomes larger than a value R δ,c1 which is approximately 85. This theoretical result provides an explanation of the laboratory observation of disturbances of the laminar regime which appear when the Reynolds number becomes larger than a value which depends on the experimental apparatus but ranges around 100.
(ii) When the Reynolds number becomes larger than R δ,c1 , but it stays smaller than approximately 200, the instability is restricted to small parts of the oscillation cycle, close to the inversions of the velocity of the oscillating plate, from t = t 1 + nπ to t = t 2 + nπ. This hydrodynamic regime can be paired with the 'disturbed laminar' regime.
(iii) When the Reynolds number is increased beyond 200 the range of unstable phases of the cycle rapidly increases until for R δ = R δ,c2 550, the instability pervades a much larger part of the cycle and in particular the whole accelerating phases. In this range of the Reynolds number, the dynamics of the perturbations is expected to give rise to the 'intermittently turbulent regime'.
(iv) Finally, when the Reynolds number becomes larger than R δ,c3 850, the instability spreads all over the cycle. Hence, turbulence is expected to be present during the whole cycle and the 'fully turbulent' regime to be observed (see the experimental measurements of Hino et al. (1983) and Jensen et al. (1989)).
In figure 11, the results of the present analysis are plotted along with the data of Hino et al. (1976), who investigated the transition process in the oscillatory boundary layer which is generated in a pipe of circular cross-section. Hence, the flow regimes observed by Hino et al. (1976) depend not only on the Reynolds number but also on the ratio λ = R * /δ * , R * being the radius of the pipe section. Even though different terms were used by Hino et al. (1976), namely (i) 'laminar flow at low Reynolds numbers', (ii) 'weakly turbulent flow at transition Reynolds numbers', (iii) 'conditionally turbulence at higher Reynolds numbers', the velocity signals plotted in figures 3, 4 and 6 of the paper of Hino et al. (1976) suggest that the flow regimes (ii) and (iii) can be associated with the disturbed laminar regime and intermittently turbulent regime, respectively. As expected, the results of the present stability analysis tend to provide a fair prediction of the flow regimes observed during the experiments of Hino et al. (1976) when λ tends to be much larger than unity, even though further experiments would be necessary to fully support the theoretical analysis.
As already pointed out, nonlinear effects certainly affects the dynamics of the perturbations when their amplitude becomes large, and quantitative differences are expected to be present when the present results are compared with the experimental observations. Indeed, differences amongst the critical values of the Reynolds number predicted by the present stability analysis and those observed in the laboratory experiments are present. Nevertheless, the present results provide a physical insight into the mechanisms giving rise to the different flow regimes characterizing an oscillatory boundary layer. Only direct numerical simulations (DNS) can be used to follow the growth of the flow perturbations until they attain their mature stage (see, amongst others, Akhavan, Kamm & Shapiro (1991b), Verzicco & Vittori (1996), Vittori & Verzicco (1998), Costamagna, Vittori & Blondeaux (2003) and Mazzuoli, Vittori & Blondeaux (2011)). In particular, Thomas et al. (2014) and Ramage et al. (2020) made DNS of oscillatory Stokes layers and showed that perturbations of the laminar regime grow rapidly after the wall shear stress reverses its direction. However, because of the high computational costs, DNS do not allow a detailed investigation of the parameter space which would require a large number of simulations for different values of the Reynolds number (R δ ) and of the size of the computational domain (α). Hence, the results of this linear stability analysis could help to select the numerical simulations to be made.