aa r X i v : . [ phy s i c s . f l u - dyn ] F e b Nonmodal Tollmien-Schlichting waves
JorisC.G.Verschaeve ∗ Noen vei 98, 0765 Oslo, NorwayFebruary 25, 2020
Abstract
The instability of flows via two-dimensional perturbations is analyzedtheoretically and numerically in a nonmodal framework. The analysis isbased on results obtained in [Verschaeve et al. (2018)] showing the invis-cid character of the growth mechanism of these waves. In particular, it isshown that the formulation of this growth mechanism naturally reducesto the eigenvalue problem for the energy bound formulated by [Davis andvon Kerczek (1973)]. The eigenvalue equation by [Davis and von Kerczek(1973)] thus allows for a broader interpretation. It provides the discretegrowth rates for the base flow in question. In addition to this eigenvalueproblem, a corresponding eigenvalue problem for the phase speed of theperturbations can be extracted from the equations found in [Verschaeve et al. (2018)]. These two eigenvalue equations relate to the Hermitianand skew-Hermitian part, respectively, of the nonmodal equations, cf.[Schmid (2007)]. In contrast to traditional Orr-Sommerfeld modal analy-sis, the above eigenvalue equations define an orthogonal set of eigenfunc-tions allowing to decompose the perturbation into base perturbations withdiscrete growth rates and frequencies. As a result of this decomposition,it can be shown that the evolution of two-dimensional perturbations isgoverned by two mechanisms: A first one, responsible for extracting andreturning energy from and to the base flow, in addition to viscous dis-sipation and, a second one, responsible for dispersing energy among thedifferent base perturbations constituting the perturbation. As a generalresult, we show that the stability of a flow is not only determined by thegrowth rates of the base perturbations, but it is also closely related to itsability to disperse energy away from the base perturbations with positivegrowth rates to the ones with negative growth rates. We illustrate theabove results by means of three well known shear flows, Couette flow,Poiseuille flow, and the boundary layer flow under a solitary wave.
Focus on two or three dimensional perturbations has shifted throughout thehistory of research on hydrodynamic stability. Whereas, as a consequence of ∗ Email address for correspondence: [email protected]
This section presents the basic equations on which the present analysis in sec-tion 3 is founded.In the present treatise we shall consider steady and unsteady base flows U base in horizontal direction, with the wall normal direction in z : U base = U ( z, t ) e x . (1)We introduce a perturbation velocity u ′ = ( u ′ , v ′ , w ′ ) in the streamwise, span-wise and wall normal direction, defined by: u ′ = ( u ′ , v ′ , w ′ ) = ( u ns , v ns , w ns ) − ( U ( z, t ) , , , (2)where ( u ns , v ns , w ns ) satisfies the Navier-Stokes equations. The energy of theperturbation is given by: E p = 12 Z V u ′ + v ′ + w ′ dV, (3)which is integrated over the entire volume of interest V . For time dependentflows, [6] derived a bound for the perturbation energy for the Navier-Stokesequations: E p ( t ) E p ( t ) ≤ exp t Z t µ ( t ′ ) dt ′ , (4)where µ is the largest eigenvalue of the following linear system:1 Re ∆ u ′ − S base ( t ) · u ′ − ∇ p = 12 µ u ′ (5) ∇ · u ′ = 0 , (6)3here the tensor S base is the rate of strain tensor given by the base flow, equation(1). Compared to the equations in [6], a missing factor of 1 / t , the eigenvalue µ isalso time dependent. If µ < x and y , we consider a single Fourier component of u ′ : ( u ′ , v ′ , w ′ )( x, y, z, t ) = ( u, v, w )( z, t ) exp i ( αx + βy ) . (7)This allows us to eliminate p from the equations (5-6), resulting into1 Re L w + i α (cid:8) D U w + 2
DU Dw (cid:9) + i β DU ζ = 12 µ L w, (8) − Re L ζ − i β DU w = 12 µ ( − ζ ) (9)where L is the Laplacian defined by: L = D − k . (10)We introduced the following shorthand notations: D = ∂∂z , (11) k = α + β . (12)The system of four equations (5-6), has been reduced to two, by means of thenormal vorticity component ζ : ζ = i ( αv − βu ) . (13)As mentioned above, the system (5)-(6) and the system (8)-(9) result from abound on the energy of the perturbation starting from the Navier-Stokes equa-tions. As we shall see in section 3, an equivalent eigenvalue problem can beobtained by introducing a parabolized stability equation ansatz for the pertur-bation and searching for the maximum growth rate.Nonmodal stability analysis is based on the linearized Navier-Stokes equa-tions, which can be written in the present setting as follows, (cid:18) ∂∂t + i αU − Re L (cid:19) L w − i αwD U = 0 , (14) (cid:18) ∂∂t + i αU − Re L (cid:19) ζ − i βwDU = 0 . (15)We refer to [14, 13] for a thorough derivation of equations (14) and (15). Givenan initial perturbation ( w , ζ ) at time t , equations (14) and (15) can be inte-grated to obtain the temporal evolution of ( w, ζ ) for t > t . Nonmodal theory4ormulates the stability problem as finding the initial condition ( w , ζ ) maxi-mizing the perturbation energy E p ( t ) of ( w, ζ ) at time t > t . This perturbationenergy E p is the sum of two contributions, one from the wall normal component w and one from the normal vorticity component ζ : E p ( t ) = E w ( t ) + E ζ ( t ) = 12 b Z a k | Dw | + | w | dz + 12 b Z a k | ζ | dz, (16)where a = − b = 1 for the enclosed base flows and a = 0 and b = ∞ forthe boundary layer flow in the present discussion.The optimization problem can then be formulated by maximizing E p for aperturbation ( w, ζ ) satisfying (14) and (15) and having an initial energy E p ( t ).One way of solving this optimization problem is by means of the adjoint equationas in [11]. Another approach for finding the optimal perturbation, which isemployed in the present treatise, consists in formulating the discrete problemfirst and computing the fundamental solution matrix X ( t, t ) of the system ofODEs, cf. references [18, 14, 13] for details. The energy E p is then given interms of X and the initial condition. Details of the implementation are given in[20, appendix A]. By computing E p ( t ) one way or the other, we can compute theamplification G from time t to t of the optimal perturbation for wave numbers α and β : G ( α, β, t , t, Re ) = max ( w ,ζ ) E p ( t ) E p ( t ) . (17)We remark that the initial condition ( w , ζ ) from which the optimal pertur-bation starts, might be different for each point in time t , when tracing G as afunction of t , cf. section 3. The maximum amplification G max ( Re ), which canbe reached for a given Reynolds number Re , is obtained by maximizing G overtime, initial time and wavenumbers: G max = max α,β,t ,t G. (18)In the following, we shall distinguish between the following three types of per-turbations: • Streamwise streaks.These are perturbations independent of the streamwise coordinate x . Theycan be computed by setting α = 0. For this case, the normal velocitycomponent w is a slowly decaying function entering equation (15) as asource term together with U . • Two dimensional perturbations.If β = 0, the governing equations (14) and (15) are decoupled. All possiblegrowth is restricted to w , as the normal vorticity ζ only displays decay. • Oblique perturbations.These are all perturbations with α = 0 and β = 0.5he present treatise focuses on perturbations of the second type which we shallcall nonmodal Tollmien-Schlichting waves due to their wave-type character. Thethird perturbation type can be considered superpositions of nonmodal Tollmien-Schlichting waves and streamwise streaks. For this type of perturbation, thenonmodal Tollmien-Schlichting wave travels in the plane spanned by the wavenumber vector k = ( α, β ) and e z .For streamwise streaks ( α = 0), on the other hand, a scaling argument as in[8, 14] allows us to rewrite equations (14) and (15): (cid:18) ∂∂τ − L (cid:19) L w = 0 , (19) (cid:18) ∂∂τ − L (cid:19) ˜ ζ − i βw ∂∂z U = 0 , (20)where τ = t/ Re is a slowly varying time scale. The normal vorticity ˜ ζ is scaledby Re : ˜ ζ = 1 Re ζ ( z, τ ) . (21)Equation (19) corresponds to slow viscous damping of w , which also holds forthe homogeneous part of equation (20) for ˜ ζ . The second term in (20) representsa forcing term which varies on the temporal scale of the base flow. Therefore,streamwise streaks display temporal variations on the time scale of the base flow.A similar scaling argument, however, does not hold for nonmodal Tollmien-Schlichting waves ( α > w is decomposed into ashape function ˜ w and an exponential factor: w = ˜ w ( z, t ) exp − i t Z t Ω( t ′ ) dt ′ , (22)where the real part of Ω accounts for the oscillatory character of w and theimaginary part of Ω is the growth rate of the perturbation. We remark that in[20], the factor − i has been absorbed into the definition of Ω.In order to define the shape function ˜ w univocally, all growth is restricted toΩ. Somewhat different to [1], a normalization condition is defined on the entire6inetic energy ˜ E w of the shape function ˜ w :˜ E w = 12 b Z a k | D ˜ w | + | ˜ w | dz. (23)The normalization constraint can thus be written as: b Z a ˜ w † L ∂ ˜ w∂t dz = 0 , (24)from which it follows that we can require that ˜ E w is unity for all times˜ E w = − k b Z a ˜ w † L ˜ w dz = 1 . (25)Equation (14) becomes then: ∂ t L ˜ w − iΩ L ˜ w = 1 Re L ˜ w + i α (cid:0) D U − U L (cid:1) ˜ w (26)Multiplying by ˜ w † and integrating in z , leads to a formula for Ω:Ω = − i2 k Re b Z a ˜ w † L ˜ w dz + α k b Z a ˜ w † D U ˜ w − ˜ w † U L ˜ w dz (27)In order to facilitate the notation, we shall write for the real and imaginaryparts of Ω: Ω = ω + i σ. (28)From equation (27), [20] derived a formula for the growth rate σ : σ = − k Re b Z a |L ˜ w | dz − i α k b Z a DU (cid:8) ˜ w † D ˜ w − D ˜ w † ˜ w (cid:9) dz (29)= − k Re b Z a |L ˜ w | dz + α k b Z a DU { ˜ w r D ˜ w i − ˜ w i D ˜ w r } dz (30)As [20] noted, the first term on the right hand side represents viscous dissipationand is always negative. The second term, however, can, depending on U and˜ w , be positive or negative. Only when this term is positive and in magnitudelarger than the viscous dissipation, growth of E w can be observed. As is im-mediately evident, the second term is multiplied by α = p k − β , which ismaximum for β = 0 for a given k . Therefore, nonmodal Tollmien-Schlichtingwaves will display larger growth rates when aligned with the base flow for cases7here the energy contribution by the normal vorticity component is negligible.This result for nonmodal Tollmien-Schlichting waves corresponds to Squire’stheorem for modal perturbations. We remark that for many flow situations, thenormal vorticity ζ does not vanish and that larger energy growth can often beobtained for E ζ than for E w . Under these circumstances the superposition (orthe pure streamwise streaks) will have a larger amplification than the nonmodalTollmien-Schlichting wave aligned with the flow.Performing some algebraic manipulations, [20] expressed σ as a function ofthe rate of strain S base for two-dimensional perturbations ( β = 0): σ = − α Re b Z a |L ˜ w | dz − b Z a (cid:0) ˜ u † , ˜ w † (cid:1) S base (cid:18) ˜ u ˜ w (cid:19) , (31)where ˜ u is the horizontal velocity component given by:i α ˜ u = − D ˜ w. (32)In the oblique case, a similar relation holds for the projection of S base onto thewave number vector k = ( α, β ), cf. [20]. Formula (31) corresponds to a tilt ofthe velocity vector (˜ u, ˜ w ) T by the rate of strain tensor S base which can be seenas the Orr-mechanism in a nonmodal framework. As [20] concluded, the growthmechanism itself is always inviscid, holding also for modal Tollmien-Schlichtingwaves which are commonly thought of as slow viscous instabilities, cf. for ex-ample [10] and [2]. Whether growth of two-dimensional perturbations is fast orslow is, as formula (31) suggests, primarily a property of the base flow profile U .In the following, we shall further develop the formalism developed in [20]and show how it allows us to deepen the theoretical understanding of nonmodalTollmien-Schlichting waves.Concerning the generation of numerical results, we employed the same nu-merical framework as in [20]. For the solution of the eigenvalue equations, adiscretization based on Shen-Legendre polynomials is employed [15]. The timedependent equations are discretized by means of Shen-Chebyshev polynomials[16] and integrated via a Runge-Kutta integrator, see [20] for implementationdetails and for verification and validation. We shall first present theoretical derivations in section 3.1 before going over tonumerical results in section 3.2. 8 .1 Theoretical considerations
In the following, we shall often use the energy scalar product for two functions φ and ψ satisfying the boundary conditions at a and b : h φ, ψ i = 12 k b Z a Dφ † Dψ + k φ † ψ dz. (33)Based on the definition of ˜ w , equation (22), the energy E w of the perturbationcan be expressed by means of the growth rate σ : E w = 12 b Z a k | D ˜ w | + | ˜ w | dz exp 2 t Z t σ ( t ′ ) dt ′ = exp 2 t Z t σ ( t ′ ) dt ′ . (34)This allows us to find a bound on E w by searching for the shape function ˜ w which maximizes σ for any point in time t , ie. we have the following variationalproblem:max ˜ w σ = max ˜ w − k Re b Z a L ˜ w † L ˜ w dz + α k b Z a DU { ˜ w r D ˜ w i − ˜ w i D ˜ w r } dz. (35)However, ˜ w has to satisfy the constraint that its energy is unity, which leads tothe following Lagrangian: L γ [ ˜ w r , ˜ w i , λ ] = − k Re b Z a L ˜ w † L ˜ w dz + α k b Z a DU { ˜ w r D ˜ w i − ˜ w i D ˜ w r } dz + λ − b Z a k | D ˜ w | + | ˜ w | dz (36)Stationarity with respect to ˜ w leads to the following eigenvalue equation:1 Re L ˜ w + i α (cid:0) DU D ˜ w + D U ˜ w (cid:1) = λ L ˜ w, (37)where the eigenvalue λ corresponds to σ , which can be seen by multiplyingequation (37) by ˜ w † and integrating in wall normal direction. It is straightfor-ward to verify that the differential operator on the left hand side of (37) is aHermitian operator. The eigenvalues are thus real and the eigenfunctions areorthogonal with respect to the energy scalar product. Seen the definition of theenergy (34), equation (37) is equivalent to the eigenvalue system (8-9), in thecase β = 0, previously found by [6]. Formulating a bound on the energy of theperturbation or alternatively formulating a bound on its growth rate via thepresent parabolized stability equation formulation produces equivalent systems,9hich supports the approach developed in [20].For Re → ∞ , the variational problem (35) reduces tomax ˜ w σ = max ˜ w α k b Z a DU { ˜ w r D ˜ w i − ˜ w i D ˜ w r } dz, (38)which, when using Euler’s formula:˜ w ( z ) = r ( z ) exp θ ( z ) , (39)becomes max r,θ σ = α k b Z a r DθDU dz. (40)At extremum, while respecting constraint (25), we obtain the following eigen-value problem: Dθ = α λ DU (41) λ (cid:0) D − k (cid:1) r + α DU ) r = 0 . (42)Using equation (41), the maximum for the growth rate can thus be written as:max r,θ σ = max r σ = max r α k vuuut b Z a r ( DU ) dz. (43)As can be seen from equation (42), in the inviscid case, the eigenvalues λ comein pairs of positive and negative values, corresponding to growth and decay re-spectively. More interestingly, however, is the fact, that the phase change Dθ ofthe eigenfunctions of system (41-42) is proportional to the rate of strain of thebase flow, cf. equation (41). This can be interpreted as a resonance mechanismof the system selecting the perturbation whose phase change matches the strainrate best.However, as we shall see in the following, equation (37) and equations (41)and (42) only describe one aspect of the physical mechanism of growth of non-modal Tollmien-Schlichting waves. In order to obtain a complete picture of thegrowth mechanism, we return to the constraint, equation (24), on which thepresent parabolized stability equation formalism is founded. Writing out the10eal and imaginary part of constraint (24), we find:12 ddt b Z a ( D ˜ w r ) + ( D ˜ w i ) + k ( ˜ w r + ˜ w i ) dz +i b Z a (cid:18) − ˜ w i L ∂∂t ˜ w r + ˜ w r L ∂∂t ˜ w i (cid:19) dz = 0 . (44)The real part of equation (44) corresponds to the conservation of energy imposedon ˜ w and which ultimately defines the growth rate σ , equation (30). On theother hand, using equation (26), the constraint on the imaginary part can bereformulated as: ω = α k b Z a
12 ( D U ) | ˜ w | + U (cid:0) | D ˜ w | + k | ˜ w | (cid:1) dz, (45)which coincides with the imaginary part of equation (27). This implies that thepresent parabolized stability equation formalism could have been derived byonly formulating a constraint on the real part of equation (44). In other words,imposing the constraint (44) on the imaginary part is redundant. Equation (45)is a formula for the frequency of nonmodal Tollmien-Schlichting waves. If thedependence of ˜ w on α had been known, equation (45) would give us an explicitdispersion relation for the frequency ω as a function of the wavenumber α .Given equation (45), the phase speed c of the nonmodal Tollmien-Schlichtingwave can be written as: c = ωα = 12 k b Z a
12 ( D U ) | ˜ w | + U (cid:0) | D ˜ w | + k | ˜ w | (cid:1) dz. (46)It consists of two terms. The first one can be seen as a weighted integral of D U , with | ˜ w | / k being the weight function. The second term is a weightedintegral of U with the energy density of ˜ w as a weight. As the energy densityof ˜ w sums to unity, the weighted integral of U gives in fact a mean velocity.This term adds more weight to the value of U where the energy density of theperturbation is largest. As such the second term indicates that the nonmodalTollmien-Schlichting wave travels with a wave speed comparable to the baseflow velocity, if the influence of the first term is negligible.Equation (45) can also be interpreted as the solution to the constrainedoptimization problem defined by the following Lagrangian: L [ ˜ w r , ˜ w i , ν ] = α k b Z a
12 ( D U ) | ˜ w | + U (cid:0) | D ˜ w | + k | ˜ w | (cid:1) dz + ν − b Z a k | D ˜ w | + | ˜ w | dz . (47)11e remark that ν is not the kinematic viscosity, but a Lagrangian multiplierin order to satisfy constraint (25). At stationarity, we obtain the followingeigenvalue equation: α (cid:18) D U ˜ w − DU D ˜ w − U L ˜ w (cid:19) = − ν L ˜ w. (48)Multiplying this equation by ˜ w † and integrating in z , we find that at extremum ν corresponds to ω . A straightforward calculation allows to verify that equation(48) is self-adjoint and therefore its eigenvalues are real and its eigenfunctionsorthogonal with respect to the energy scalar product.Turning again to Euler’s formula, equation (39), we obtain for the eigenvalueproblem (48): Dθ = 0 (49) α (cid:0) D U r − D ( U Dr ) + 2 k U r (cid:1) = ν (cid:0) − D r + k r (cid:1) (50)Equation (49) corresponds to the fact, that the system is self-adjoint (in thesense that it can be written as a real and not a complex equation), opposedto the above truly Hermitian one, equation (37). Therefore, contrary to theeigenfunctions of equation (37), the eigenfunctions of (48) can be normalized tobe real, ie. θ = 0. Equation (49) also contrasts with the formula for the phaseof the eigenfunctions with maximum growth rate in the inviscid case, equation(41), indicating that in general the optimum of one system cannot be the opti-mum of the other.Eigenvalue problem (37) is related to the buckling problem of a thin rod[19]. The eigenvalues λ i , equation (37), appear as a decreasing sequence: λ > λ > . . . > λ i > . . . . (51)As mentioned in [6], growth of nonmodal Tollmien-Schlichting waves is possibleonly if at least λ is larger than zero. As can be seen from equation (35), theviscous contribution, ie. the first term on the left hand side of (35), cannotbe bounded from below, leading to discrete eigenvalues towards −∞ , at leastfor bounded domains. An upper bound λ max > λ is given in [6, formula(2.8)] for bounded domains. However, even for unbounded domains, we canuse equation (43) to find an upper bound for the second term in (35), as longas R ba ( DU ) dz < ∞ , which is the case for all flows considered in the present12iscussion. Using H¨olders inequality, we find: λ ≤ α k b Z a r ( DU ) dz (52) ≤ α k b Z a r dz b Z a ( DU ) dz (53) ≤ α k b Z a ( DU ) dz | {z } =: λ , (54)where we have used the fact that the energy of the perturbation is unity andwe thus have the following bound for the integral of r :12 b Z a r dz = 12 b Z a | ˜ w | dz ≤ k b Z a | D ˜ w | + k | ˜ w | dz = 1 . (55)In general, the mathematical theory on higher order Sturm-Liouville problemswith a second order differential operator in the term multiplied by the eigenvalueis rather limited [5]. A mathematical proof of property (51) of the discretespectrum of equation (37) is still an open question. In addition, for unboundeddomains, a mathematical analysis of its discrete and continuous spectrum is stillmissing. However, these questions surpass the scope of the present discussion.The numerical results in 3.2 give an indication that a discrete spectrum withproperty (51) exists under some conditions even for unbounded domains.Equation (48) in this respect is even more exceptional in the sense thatthe term multiplied by the eigenvalue is as the left hand side, a second orderdifferential operator. The eigenvalues of (48) are observed to lie in a specificrange, cf. section 3.2: ν i ∈ [ ω min , ω max ] . (56)From a physical point of view, this is sensible, as it should not be possible for anonmodal Tollmien-Schlichting to travel with infinite speed. Property (56) canbe proven the following way.For the boundary layer flow in the present discussion, we havelim z →∞ U ( z, t ) = U b ( t ) , (57)where U b is not necessarily zero. In this case the integral R ∞ | U | dz is not defined.However, we assume that when subtracting U b from U , its integral is bounded,ie. we have: ∞ Z | U − U b | dz < ∞ . (58)13sing equation (45), we can find a bound for | ω | :2 k α | ω | ≤ (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) b Z a
12 ( D U ) | ˜ w | + U (cid:0) | D ˜ w | + k | ˜ w | (cid:1) dz (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) (59) ≤ (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) b Z a
12 ( D U ) | ˜ w | dz (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) + (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) b Z a U (cid:0) | D ˜ w | + k | ˜ w | (cid:1) dz (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) (60)= (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) b Z a
12 ( D U ) | ˜ w | dz (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) + (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) b Z a ( U − U b + U b ) (cid:0) | D ˜ w | + k | ˜ w | (cid:1) dz (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) (61) ≤ b Z a (cid:12)(cid:12) D U (cid:12)(cid:12) dz b Z a | ˜ w | dz + b Z a | U − U b | dz b Z a (cid:0) | D ˜ w | + k | ˜ w | (cid:1) dz + | U b | b Z a (cid:0) | D ˜ w | + k | ˜ w | (cid:1) dz (62) ≤ b Z a (cid:12)(cid:12) D U (cid:12)(cid:12) dz + 2 k b Z a | U − U b | dz + 2 k | U b | , (63)where we have made use of H¨older’s inequality. A consequence of the bound-edness of the eigenvalues of (48) is that we do not have a hierarchy of theeigenvalues as in (51). This manifests itself during the numerical calculations inthe fact that when increasing the numerical resolution, the eigenvalues for finerresolutions appear in between the ones for coarser resolutions.On the other hand, as the eigenvalues of equation (37) appear in a decreas-ing sequence, equation (37) allows us to define a proper set of orthonormaleigenfunctions with respect to the energy scalar product: { λ i , φ i } , i = 0 , , , . . . . (64)In the following, we shall call the functions φ i VKD-modes in honor of Von Ker-czek and Davis who were the first to derive eigenvalue system (5-6). Prescindingfrom the continuous spectrum in the case of unbounded domains, the nonmodalTollmien-Schlichting wave can then be expanded onto the VKD-modes: w = X i =0 c i φ i , (65)where c i = h φ i , w i . (66)In quantum mechanics, the VKD-modes would correspond to the base states ofa quantum mechanical system with | c i | its probability to be found in this state.14he analogy is not complete, as in the present case | c i | is the energy of theVKD-mode. However, the VKD-modes can be seen as the base perturbations ofthe governing system (14). Defining c = ( . . . , c i , . . . ) t as the coefficient vector,we can write the governing equation (14) in Heisenberg form: ddt c + Fc = Λc − i Nc . (67)The matrices in equation (67) have the following meaning. Matrix F accountsfor the temporal change of the VKD-modes. Its elements are defined by F ij = (cid:28) φ i , ∂∂t φ j (cid:29) . (68)The base perturbations φ i result from an eigenvalue problem for each point intime. Therefore, they display temporal variations on the same scale as the baseflow. In order to trace φ i in time, a mapping of the VKD-modes between differ-ent points in time needs to be established. This poses some challenges as shallbe elaborated more in detail in appendix A. For steady base flows the matrix F vanishes. Many boundary layer flows are governed by two well separated timescales. For the boundary layer flow under a solitary wave, the base flow varieson a slow time scale, namely Re /
2, whereas the nonmodal Tollmien-Schlichtingwave varies on a fast time scale. As F varies on the same scale as the base flow,we can neglect it for higher Reynolds numbers. This is similar to neglectingnonparallel effects in boundary layers developing in streamwise direction. In or-der to ease the discussion, we shall only consider the cases where F is negligiblein the following. The governing equation becomes then ddt c = Λc − i Nc . (69)The matrix Λ is a diagonal matrix with the growth rates on its diagonal:Λ ij = λ i δ ij . (70)The elements of matrix N are defined as follows: N ij = α k b Z a D U φ † i φ j + U (cid:16) Dφ † i Dφ j + k φ † i φ j (cid:17) dz. (71)As N is a Hermitian matrix, the matrix i N is skew-Hermitian, and since Λ isHermitian, the right hand side of equation (69) corresponds to the decomposi-tion of (14) into its Hermitian and skew-Hermitian part: ∂∂t L w = 1 Re L w + i α (cid:18) DU Dw + 12 D U w (cid:19)| {z }
Hermitian + i α (cid:18) D U w − DU Dw − U L w (cid:19)| {z } skew − Hermitian . (72)15eturning to equation (69), matrix Λ can be interpreted as a base Hamiltonianof the system, whereas − i N acts as a perturbation potential. In this view thebase perturbations are growing or decaying exponential functions under actionof a conservative potential. If we had chosen the eigenfunctions of (48), N would be diagonal. In this picture i N would be the base Hamiltonian and Λ aperturbation potential. The base perturbations in this case would be neutrallystable waves. This choice of discretization and the approximation by Feynmanpath integrals of the fundamental solution shall be discussed in more detailin appendix B. For the above mentioned reasons, we shall only consider theeigenfunctions of (37) in the following. Equation (69) can be slightly rewritten: ddt c = Lc − i Mc . (73)The matrix L is still diagonal but with its elements defined by L ij = ( λ i − i ω i ) δ ij , (74)where ω i = α k b Z a D U | φ i | + U (cid:0) | Dφ i | + k | φ i | (cid:1) dz (75)is the diagonal element of N . With L as the base Hamiltonian, the base pertur-bations would be growing or decaying waves with frequency ω i given by (75).Independent on the view adopted, the growth rate σ of the resulting nonmodalTollmien-Schlichting wave is given by: σ = 12 1 E w dE w dt = 12 d c † dt c + c † d c dt c † c = c † Λcc † c . (76)Likewise the frequency ω of the nonmodal Tollmien-Schlichting wave can becomputed by: ω = c † Ncc † c . (77)The decompositions (69) or (73) suggest a first approximation by neglectingthe skew-Hermitian part, ie. the matrices N or M , respectively. As the VKD-modes are uncoupled in this case, the optimal perturbation would be identicalto the zeroth VKD-mode, displaying extreme growth for rather small Reynoldsnumbers. This contrasts with observations for many flows, such as Couetteflow, which display only weak growth for two-dimensional perturbations. Onthe other hand, flows with adverse pressure gradients, such as the boundarylayer flow under a solitary wave, have a tendency to favor two-dimensionalperturbations.Equation (69) provides a way to understand the effective growth of non-modal Tollmien-Schlichting waves. As mentioned above, in order for growth ofnonmodal Tollmien-Schlichting waves to be possible, at least the largest eigen-value λ needs to be positive. In order to simplify the following discussion, we16hall assume that only λ is positive and all other eigenvalues λ i , i = 1 , . . . arenegative. In this case, the only VKD mode able to extract energy from the baseflow is φ . The first line of equation (69) reads as follows: ddt c = λ c − i X j =0 N j c j . (78)As λ is the only positive eigenvalue, the first term on the right hand side (78)is the only production term present in the system. The second term on the righthand side is a dispersion term distributing energy from the zeroth VKD modeto other VKD modes or vice versa. Therefore, even if the numeric value of λ is large, dispersion between VKD modes can lead to a drain of energy towardsother VKD modes, which stabilizes the system, as these modes dissipate energy.As we shall see in section 3.2, the dispersion term in equation (78) is of principalimportance for the evolution of nonmodal Tollmien-Schlichting waves. A boundfor this term can be found by: | X j N j c j | = | − α b Z a k D U φ † w + 12 k U (cid:16) Dφ † Dw + k φ † w (cid:17) dz | (79) ≤ α k | b Z a | D U φ † w − D ( U Dφ † ) w + k U φ † w | dz | (80) ≤ α k b Z a | D U φ † − D ( U Dφ † ) + k U φ † | dz b Z a | w | dz (81) ≤ α k b Z a | D U φ † − D ( U Dφ † ) + k U φ † | dz × k b Z a | Dw | + k | w | dz. (82)This allows us to define a measure m , max of the dispersion of the most dangerousVDK mode for the base flow U in question: m , max = max c P i,j c † i N i N j c j P i c † i c i ≤ α k b Z a | D U φ − D ( U Dφ ) + k U φ | dz (83)The left hand side of inequality (83), ie. the measure m , max , can easily befound numerically by computing the spectral radius of the matrix with ele-ments N i N j . In principle, the analytic bound on the right hand side of (83)gives us a possibility to relate the shape of the base flow profile U to its stabilityproperties. Its interpretation is, however, not straightforward as φ implicitly17s a function of U , resulting from eigenvalue system (37). In principle, base flowprofiles U with smaller bounds would be more unstable than those with largerbounds. A crude heuristic argument would suggest that profiles with oscillatorybehavior around zero would probably be good candidates for unstable flows asa simple substitution of U by exp ± i √ kz reduces the argument of the integralto a single term. On the other hand, this heuristic argument would also implythat the perturbation with wavenumber k bringing D U/ k U as close tozero as possible would experience larger growth.Another estimate of the dispersion properties of the zeroth VDK mode isgiven by the temporal change of its energy, | c | . From equation (78), we obtain: ddt | c | = 2 λ | c | + i X j =0 c † M k c j − c † j M † j c (84)= 2 λ | c | + i X j =0 c † M j c j − c † j M j c (85)= 2 λ | c | − X j =0 imag (cid:16) c † M j c j (cid:17) (86)Besides the measure m , max , which is purely a property of the base flow U ,we shall investigate the following quantities m and n , being actual dispersionmeasures for a given nonmodal Tollmien-Schlichting wave w : m = P i,j c † i N i N j c j P i c † i c i (87)= α k E w | b Z a D U φ † w + U (cid:16) Dφ † Dw + k φ † w (cid:17) dz | (88) n = i X j =0 c † M j c j − c † j M j c (89)= i α k h φ , w i b Z a D U w † φ + U (cid:0) Dw † Dφ + k w † φ (cid:1) dz − (90) h w, φ i b Z a D U φ † w + U (cid:16) Dφ † Dw + k φ † w (cid:17) dz (91)In the following, we shall continue the present investigation by means of a nu-merical analysis of three shear flows. 18 .2 Numerical results Optimal perturbations for Couette flow have been investigated by [3]. Theyfound that for Couette flow at Re = 1000, the global maximum is reached at t =117 with G = 1184 . α = 0 .
035 and β = 1 . β = 0, we observe that optimalnonmodal Tollmien-Schlichting waves display only decay. Perturbations with β = 0, have been found by [3] to display only weak growth. The maximumamplification reached by a perturbation with β = 0, they found, is at t = 8 . G = 13 . α = 1 .
21. However, as canbe seen from figure 2, even in this case the perturbation with α = 1 . β = 2 . α = 1 .
21 and β = 0 whose initial conditionleads to the maximum at t = 8 .
7. The energy E w of this Tollmien-Schlichtingwave is plotted together with the amplification G of the optimal perturbationfor α = 1 .
21 and β = 0 in figure 3. In the following, the initial energy E w ( t )at time t is without loss taken to be unity. In figure 4, the growth rate σ ,equation (30), of this nonmodal Tollmien-Schlichting wave is plotted. As canbe observed the growth rate is bounded by λ , but comes close to this value ataround t = 5 .
5. At approximately the same time the dispersion measure n ,equation (88) changes sign indicating that energy is first transferred from otherVDK modes to the zeroth mode and then back again. This coincides well withthe energy fraction contained in the zeroth VKD mode, see bottom figure 4.First n is positive, allowing to accumulate energy in the zeroth VKD modeand to extract more energy from the base flow. However, as the energy fractionin the zeroth VKD mode grows, n diminishes before crossing sign and invertingthe energy transfer. In figure 5, the first four eigenvalues λ i are plotted as afunction of α . As mentioned above, in the present discussion the eigenvalues λ i are enumerated in decreasing order. This implies that when plotting theeigenvalues as a function of, for example, the wave number α , the curves candisplay kinks as for the eigenvalues λ and λ at α = 0 .
18, where a modecontinuation approach would have suggested that the growth rate of one modesurpasses the one of another.For α = 1 .
21, we observe that the first three VKD modes have a positive19rowth rate λ i . This fits the above picture, cf. figure 4 bottom, where energyis transferred mainly between the first three growing VKD modes in order toobtain an optimal amplification at t = 8 .
7. If, we instead consider the optimalperturbation with maximum at t = 8 . α = 0 .
3, which, cf. figure 5, hasonly a single positive eigenvalue λ , we observe a somewhat different picture.The energy is mainly concentrated in the zeroth mode.In both cases, we observe that the dispersion measure m stays well belowits bound m , max .Reaching optimal amplification at some point in time means thus to optimizethe transfer of energy between VDK modes in order to maximize the energyfraction contained in the growing VKD modes.In figure 7, the real and imaginary part of the first three VKD modes with α = 1 .
21 and β = 0 for Couette flow at Re = 1000 are plotted. The VKD modesare determined up to a constant phase, which in the present case is chosen suchthat the real part is symmetric around the origin, whereas the imaginary partis antisymmetric.Thanks to its simple form, we can compute some of the theoretical quanti-ties analytically for Couette flow. The upper bound λ max on the growth rate,equation (54), is for β = 0 equal to unity: λ max = 1 . (92)Couette flow also allows to solve the eigenvalue equation (42) for inviscid growthanalytically. In fact, eigenvalue equation (37) for viscous growth is also solvableanalytically. However, due to the complexity of the solution of a quartic equa-tion, the analytic solution becomes cumbersome. The solution of (42), can bewritten as follows: λ n, ± = α α + γ n ) , (93) γ n = π n ) , (94) r n = (cid:26) A n cos γ n z for n even A n sin γ n z for n odd (95)We obtain thus for the zeroth VKD mode: λ = α √ α + π (96) φ = A e i α λ z cos π z (97)Up to a phase, the constant A is determined by the constraint that the energyof φ is unity: | A | = 4 α α + π = 4 λ . (98)For large values of α , the inviscid growth rate λ converges towards 1 /
2. Onthe other hand, as can be observed from figure 8, the growth rate λ becomes20egative for large values of α in the viscous case. The larger the Reynoldsnumber, the longer λ stays positive. This is in accordance with formula (30),where the viscous part scales as α for large values of α . This dissipation ofsmall scales is absent in the inviscid case.In the inviscid case, we can compute the inviscid dispersion bound analyti-cally with β = 0 for Couette flow:12 α b Z a | D U base φ − D ( U base Dφ ) + k U base φ | dz = 24 α + π (cid:18) π (cid:18) π (cid:19) + α (cid:0) π − (cid:1) + 4 α (cid:18) − π (cid:19)(cid:19) (99)From formula (99), we observe that for α → ∞ , the bound displays a quadraticbehavior with α . This behavior can also be seen for m , max in the viscouscase, cf. figure 5. As such λ is for most Reynolds numbers growing fasterfor smaller α than m , max , before m , max overtakes λ , cf figure 5, indicatingthat for a certain value of α , growth and dispersion to other VKD modes reachbreak even. It is thus rather the quadratic growth of m , max than the viscousdissipation of small scales, which is at the origin of a finite value of α for whichthe optimal perturbation reaches a maximum amplification, for example thevalue α = 1 .
21 for Couette flow at Re = 1000.21 − − − − − G . . . . . . β . . . . . . α . . . . . G Figure 1: Isolines of G for the optimal perturbation for Couette flow at Re =1000 at t = 117. The plots to the left and below the contour plot show a slicealong the β - and α -axes, respectively. 22 . . . . . . . . . . . . . G β α G Figure 2: Isolines of G for the optimal perturbation for Couette flow at Re =1000 at t = 8 .
7. The plots to the left and below the contour plot show a slicealong the β - and α -axes, respectively. 23 GE w Figure 3: Amplification G of the optimal perturbation and temporal evolutionof the energy E w of the perturbation leading to a maximum amplification at t = 8 . α = 1 .
21 and β = 0 for Couette flow at Re = 1000.24 λ σ n /E w m m
0, max c /E w c /E w c /E w c /E w Figure 4: Temporal evolution of characteristic quantities of the nonmodalTollmien-Schlichting wave with α = 1 . β = 0, for Couette flow at Re = 1000.Top: Growth rate σ with upper bound λ and dispersion measures m , n and m , max . Bottom: Energy contained in the first four VKD modes.25 .0 0.5 1.0 1.5 2.0α−0.20.00.20.40.60.81.0 λ λ λ λ m
0, max
Figure 5: Growth rates λ i , equation (51), for Couette flow at Re = 1000 anddispersion measure m , max , equation (83).26 λ σ n /E w m m
0, max c /E w c /E w c /E w c /E w Figure 6: Temporal evolution of characteristic quantities of the nonmodalTollmien-Schlichting wave with α = 0 . β = 0, for Couette flow at Re = 1000.Top: Growth rate σ with upper bound λ and dispersion measures m , n and m , max . Bottom: Energy contained in the first four VKD modes. −0.5 −0.4 −0.3 −0.2 −0.1 0.0 0.1 0.2 0.3−1.0−0.50.00.51.0 z real φ real φ real φ (a) Real part −0.6 −0.4 −0.2 0.0 0.2 0.4 0.6−1.0−0.50.00.51.0 z imag φ imag φ imag φ (b) Imaginary part Figure 7: Real and imaginary part of the first three VDK modes, with α = 1 . β = 0 for Couette flow at Re = 1000.27 α −0.2−0.10.00.10.20.30.40.5 λ Re = ∞ Re = 5000Re = 1000Re = 500Re = 100
Figure 8: Growth rate λ as a function of α for Couette flow at different Reynoldsnumbers. 28 .3 Poiseuille flow Similar to Couette flow, [3] found that streamwise streaks are the dominantperturbations for Poiseuille flow. At a Reynolds number of Re = 5000, [3]calculated that the streamwise streak with α = 0 and β = 2 .
044 reaches theglobal maximum at t = 379, with G = 4897. In figure 9, the amplification G of the optimal perturbation at t = 379 is plotted in Fourier space for Poiseuilleflow at Re = 5000. As can be observed from this plot, streamwise streaksdominate over two-dimensional perturbations which display just a small peakaround α ≈ .
05. According to [3], the global maximum for a two-dimensionalperturbation ( β = 0), is reached for the nonmodal Tollmien-Schlichting wavewith α = 1 .
48 at time t = 14 .
1. As for Couette flow, this maximum is only asaddle point in Fourier space, cf. figure 10. The maximum at this time is reachedby the superposition with α ≈ . β ≈
4. As before, we concentrate onthe evolution of the nonmodal Tollmien-Schlichting wave with α = 1 .
48 and β = 0, being the optimal two-dimensional perturbation at t = 14 .
1, cf. figure11. When plotting the evolution of n , cf. figure 12, we observe, as before thatfirst energy is transferred to the zeroth VKD mode before being returned toother modes again. There is, however, a difference to Couette flow. The energycontained in VKD modes with an odd index is zero (not shown). When plottingthe real and imaginary part of the first VKD modes, cf. figure 13, we observethat VKD modes with even indices are even functions and VKD modes withodd indices are odd functions. From equation (71), we can readily infer, that inthis case, the matrix N is banded, since we have for the elements of N : N ij = 0 if i + j odd . (100)For this reason, only the even indexed VKD modes are part of the optimalperturbation, as no energy transfer between odd and even indexed VKD modes ispossible. As for Couette flow, we observe a quadratic behavior of the dispersionmeasure m , max , cf. figure 14. In addition, as we can infer from figure 14, thereare three VKD-modes with positive growth rates mainly transferring energybetween each other as depicted in figure 12. When choosing the value α = 0 . λ . As visible from figure 15, we observe, similarto Couette flow, that most of the energy is contained in the zeroth mode.For the above Couette flow, the phase speed of the nonmodal Tollmien-Schlichting wave and the phase speed of the first VKD modes is zero (not shown).This is, however, not the case for Poiseuille flow at Re = 5000, cf. figure 16.The VKD modes travel at different phase speeds, equation (75). The nonmodalTollmien-Schlichting wave propagates with a phase speed close to the averagevelocity of the base flow. 29 . . . . . . . . . . G β α G Figure 9: Isolines of G for the optimal perturbation for Poiseuille flow at Re =5000 at t = 379. The plots to the left and below the contour plot show a slicealong the β - and α -axes, respectively. 30 G β α G Figure 10: Isolines of G for the optimal perturbation for Poiseuille flow at Re = 5000 at t = 14 .
1. The plots to the left and below the contour plot show aslice along the β - and α -axes, respectively.31 GE w Figure 11: Amplification G of the optimal perturbation and temporal evolutionof the amplification E w of the perturbation leading to a maximum amplificationat t = 14 . α = 1 .
48 and β = 0 for Poiseuille flow at Re = 5000.32 λ σ n /E w m m
0, max c /E w c /E w c /E w c /E w Figure 12: Temporal evolution of characteristic quantities of the nonmodalTollmien-Schlichting wave with α = 1 . β = 0, for Poiseuille flow at Re = 5000. Top: Growth rate σ with upper bound λ and dispersion mea-sures m , n and m , max . Bottom: Energy contained in even VKD modes. −0.5 −0.4 −0.3 −0.2 −0.1 0.0 0.1 0.2 0.3−1.0−0.50.00.51.0 z real φ real φ real φ (a) Real part −0.6 −0.4 −0.2 0.0 0.2 0.4 0.6−1.0−0.50.00.51.0 z imag φ imag φ imag φ (b) Imaginary part Figure 13: Real and imaginary part of the first three VKD modes with α = 1 . β = 0 for Poiseuille flow at Re = 5000.33 .0 0.5 1.0 1.5 2.0α−0.20.00.20.40.60.81.0 λ λ λ λ λ m
0, max
Figure 14: Growth rates λ i , equation (51), for Poiseuille flow at Re = 5000 anddispersion measure m , max , equation (83).34 λ σ n /E w m m
0, max c /E w c /E w c /E w c /E w Figure 15: Temporal evolution of characteristic quantities of the nonmodalTollmien-Schlichting wave with α = 0 . β = 0, for Poiseuille flow at Re = 5000. Top: Growth rate σ with upper bound λ and dispersion mea-sures m , n and m , max . Bottom: Energy contained in even VKD modes.35 t cω /αω /αω /αω /α Figure 16: Phase speed c of the nonmodal Tollmien-Schlichting wave with α =1 .
48 and β = 0 and phase speed of the four first VDK modes, for Poiseuille flowat Re = 5000. 36 .4 Boundary layer under a solitary wave Compared to the two preceding flow examples, the boundary layer flow under asolitary wave has significant distinct features. It is not only time dependent anddefined in a semi-infinite interval, but most importantly, it develops an adversepressure gradient for time t >
0. The flow is an idealization of the boundarylayer flow under a solitary wave and was originally proposed as a model flow in[17]. The outer flow of this boundary layer is given by: U outer (2 t/ Re ) = sech (2 t/ Re ) . (101)Opposed to the present discussion, we remark that in [20], time t was measuredin the scale of the outer flow. As in [20], the boundary layer equations aresolved numerically to obtain the base flow U (2 t/ Re , z ). In figure 17, the outerflow U outer and some velocity profiles for selected times are plotted. For times t < t > t > Re = 316, which, among other cases, has also been investigated in [20]. In figure18, contour plots of the amplification G in Fourier space for different times t and t are plotted. During acceleration from 2 t / Re = − t/ Re = 0, figure18(a), the only optimal perturbations displaying growth are streamwise streakswith a maximum G ≈
15. For β = 0, no growth is visible. However, whenintegrating the system (14-15) for the same duration, but in the decelerationregion from 2 t / Re = 0 to 2 t / Re = 1, figure 18(b), we observe weak growth foroptimal perturbations with β = 0 ( G ≈ . G ≈
80. However, when increasing the integration interval to 2 t/ Re = 2 and2 t/ Re = 3, figures 18(c) and 18(d), respectively, streamwise streaks show slowdecay from their peak at 2 t/ Re = 1, whereas two-dimensional perturbationsdisplay strong growth, centered around α ≈ . α,t ,t G ( α, β = 0 , t , t , Re = 316) = 3 . · , (102)where at maximum, we have α max = 0 . , (103)2 t , max / Re = 0 . , (104)2 t , max / Re = 7 . . (105)When plotting the amplification G of the optimal perturbation with α = 0 . t / Re = 0 .
509 and the evolution of the energy E w of the nonmodalTollmien-Schlichting wave with parameters given by (103-105), cf. figure 19,37e observe a qualitatively different picture than for Couette and Poiseuille flow,figures 3 and 11, respectively. The graphs of the temporal evolution of E w andof G are lying on top of each other, except for a short initial period of time, cf.figure 20 where a zoom of figure 19 is displayed. This indicates that for mostpart of the deceleration region, the energy transfer optimization from and to thezeroth VDK mode as for Couette and Poiseuille is marginal. Instead, the non-modal Tollmien-Schlichting wave evolves as if it had been an orthogonal modewith the largest growth rate. Figure 21 lends some support to this behavior.We observe that for this Reynolds number, only the zeroth VKD mode displaysregions of growth. As shall be discussed in appendix A.2, the first VKD mode,and so its eigenvalue λ , is probably already part of the continuous spectrum.The region of growth of λ is skewed towards the deceleration region of theflow. In particular, we observe that growth starts at around 2 t/ Re ≈ −
1, butstretches much further into the deceleration region. However, probably moresignificant than the extended growth in the deceleration region is the skewnessof the behavior of the dispersion measure m , max which shows the opposite be-havior. It displays much larger values in the acceleration region than in thedeceleration region. In particular, it drops below the value of λ at 2 t/ Re ≈ / t/ Re >
1. This drop of dispersion of energy to andfrom higher VKD modes in the deceleration region of the flow, in combinationwith a significant growth rate, is probably at the origin of the behavior observedin figure 19. As a result, the mode φ behaves almost as an orthogonal modeevolving independently of the other VKD modes. However, this independence isnot complete. In figure 22, several characteristic quantities are plotted for thisnonmodal Tollmien-Schlichting wave. As to be expected, σ and m honor theirupper bounds λ and m , max , respectively. On the other hand, the energy con-tent for the zeroth VKD mode evolves around 70 % of E w , indicating that, evenin the case of reduced dispersion, the Tollmien-Schlichting wave transfers, alongits course, some energy to higher VKD which for the boundary layer flow undera solitary wave at Re = 316, are thought to lie in the continuous spectrum ofequation (37). Whether this is due to a rest dispesion by the matrix N or by thematrix F , which is small in magnitude but not zero in this case, remains an openquestion. Concerning the phase speed c of the nonmodal Tollmien-Schlichtingwave, equation (46), its graph in figure 23 shows that the nonmodal Tollmien-Schlichting wave travels approximately with the phase speed ω /α of the zerothVKD mode, equation (75). The upper and lower bounds for the phase speed, ω min and ω max , respectively, cf. equation (56), are computed by searching forthe maximum and mininum eigenvalue of N . The upper bound seems to followthe speed of the outer flow, whereas the lower bound takes into account thatthe adverse pressure gradient causes a reverse flow, allowing the perturbationsto travel in opposite direction. As the continuous spectrum typically attenuatesinside the boundary layer, cf. reference [9], it is consistent that φ travels at thespeed of the outer flow. 38 z U o u t e r / U Figure 17: Inviscid outer flow U outer and profiles of the horizontal velocity com-ponent in the boundary layer under a solitary wave moving from right to left.The profiles have been multiplied by 40. The value at z = 0 of the profiles showncorresponds to the point in time t , at which the profile has been taken. Thehorizontal velocity vanishes at z = 0 in order to satisfy the no-slip boundarycondition. 39 . . . . . . . . G β α G (a) 2 t / Re = −
1, 2 t / Re = 0 . . . . . . . . . G β α G (b) 2 t / Re = 0, 2 t / Re = 1 . . . . . . . . . G β α G (c) 2 t / Re = 0, 2 t / Re = 2 . . . . . . . . G β α G (d) 2 t Re = 0, 2 t / Re = 3 Figure 18: Isolines of the amplification G ( α, β, t , t , Re δ = 316) for the bound-ary layer flow under a solitary wave, for different values of t and t . The plotsto the left and below the contour plot show a slice along the β - and α -axes,respectively. 40 GE w Figure 19: Amplification G of the optimal perturbation and temporal evolutionof the amplification E w of the perturbation leading to a maximum amplificationat 2 t/ Re = 7 .
686 with 2 t / Re = 0 . α = 0 .
369 and β = 0 for the boundarylayer flow under a solitary wave at Re = 316.41 .5 0.6 0.7 0.8 0.9 1.0 1.1 1.22t/Re02468101214 GE w Figure 20: Zoom into figure 19. Amplification G of the optimal perturbationand temporal evolution of the amplification E w of the perturbation leading to amaximum amplification at 2 t/ Re = 7 .
686 with 2 t / Re = 0 . α = 0 .
369 and β = 0 for the boundary layer flow under a solitary wave at Re = 316.42 m
0, max λ λ Figure 21: Temporal evolution of characteristic quantities of the DVK modeswith α = 0 . β = 0, for the boundary layer flow under a solitary wave at Re = 316. Growth rates λ , λ and dispersion measure m , max .43 λ σ m m
0, max c /E w Figure 22: Temporal evolution of characteristic quantities of the nonmodalTollmien-Schlichting wave with α = 0 . β = 0, for the boundary layer flowunder a solitary wave at Re = 316. 44 c ω max /αω min /αω /αω /α Figure 23: Boundary layer flow under a solitary wave with Re = 316. Displayedis the phase speed c of the nonmodal Tollmien-Schlichting wave with parametersgiven by equation (46) and phase speeds ω /α and ω /α of the first two DVKmodes, equation (75) and maximum and minimum phase speeds ω min /α and ω max /α , equation (56). 45 Conclusions
In the present treatise, we showed that the parabolized stability equation ap-proach derived in [20] leads to a Hermitian eigenvalue equation, equation (37)for the growth rate of the perturbation. The resulting set of orthonormal eigen-functions, whose temporal continuation we called VKD modes, allowed us toformulate the governing equation in Heisenberg form. The resulting matrixequation consists of a Hermitian part responsible for growth of the perturba-tion and an skew-Hermitian part redistributing energy between VKD modes.Different quantities and bounds measuring the dispersion properties betweenVKD modes have been derived. The theoretical framework developed in thepresent treatise has been applied to three shear flows, Couette flow, Poiseuilleflow and the boundary layer flow under a solitary wave. Two different regimesof a energy transfer between the VKD modes have been observed. For Couetteand Poiseuille flow, with relatively large dispersion measure m , max , equation(83), the optimal perturbation results from balancing transfer of energy fromand to growing VKD modes in such a way that growth is largest. On the otherhand, for the adverse pressure region of the boundary layer flow, a differentregime of growth became visible. As the dispersion in this case is weak, growthis almost entirely provided by energy extraction from the base flow. As we haveseen some of the energy accumulated is dispersed to the continuous spectrum.As a matter of fact, the present work leads to further questions, summarized inthe following points: • A underlying assumption of the present treatise is that growth and disper-sion break even for some specific value of α such that the resulting optimalperturbation reaches the global maximum for two-dimensional perturba-tions. Future research might find more specific conditions when and howthis happens. • A better analytic bound for m , max on the right hand side of equation(83) might furnish us with more insight on how the shape of the base flowprofile influences dispersion of energy to and from higher VKD modes. • Future research should deal with the improvement of the mathematicalfoundations of second and fourth order Sturm-Liouville problems with theeigenvalue term being a second order differential operator. In particular,semi-infinite intervals play a major role for boundary layer flows. • As the Heisenberg formulation, equation (67), allows to highlight the anal-ogy to a quantum electrodynamical system, it might be worthwhile toconsider techniques of this field to the present problem. In particular, itmight be possible to model the dispersion effect of the higher VKD modesby a mean-field theory leading to an effective growth rate.The author would like to thank Graigory Sutherland, Pawe l Wroniszewski andFlorian Schwertfirm for their support and encouragement.46
Computing VKD-modes
The numerical computation of VKD-modes requires some care. In the following,we shall elucidate two important issues when dealing with VKD-modes.
A.1 Normalizing VKD modes in time
The eigenfunctions resulting equation (37) are only determined up to a mul-tiplicative constant. In order to find the correct scaling of the eigenfunctionsin time, we return to constraint (24). By means of a semi-discrete version ofconstraint (24) at the midpoint ( t + t ) /
2, the mode φ at time t is normalizedusing its value at t :12( t − t ) b Z a (cid:0) Dφ † ( t ) + Dφ † ( t ) (cid:1) ( Dφ ( t ) − Dφ ( t ))+ k (cid:0) φ † ( t ) + φ † ( t ) (cid:1) ( φ ( t ) − φ ( t )) dz = 0(106)The real and imaginary part of equation (106) are given by: b Z a (cid:0) | Dφ ( t ) | + k | φ ( t ) | (cid:1) − (cid:0) | Dφ ( t ) | + k | φ ( t ) | (cid:1) dz = 0(107) b Z a (cid:0) Dφ † ( t ) Dφ ( t ) + k φ † ( t ) φ ( t ) (cid:1) − (cid:0) Dφ † ( t ) Dφ ( t ) + k φ † ( t ) φ ( t ) (cid:1) dz = 0(108)Equation (107) represents the conservation of energy, whereas equation (108)can be written as A − A † = 0 , (109)implying that the imaginary part of A vanishes. We can thus write: A = A r + i A i (110)= b Z a (cid:0) Dφ † ( t ) Dφ ( t ) + k φ † ( t ) φ ( t ) (cid:1) dz (111)= b Z a Dφ r ( t ) Dφ r ( t ) + Dφ i ( t ) Dφ r ( t ) + k ( φ r ( t ) φ r ( t ) + φ i ( t ) φ i ( t )) dz +i b Z a Dφ r ( t ) Dφ i ( t ) − Dφ i ( t ) Dφ r ( t ) + k ( φ r ( t ) φ i ( t ) − φ i ( t ) φ r ( t )) dz (112)After having traced the eigenfunction ˜ φ at t corresponding to the mode at t ,we can pose: φ ( t ) = e i δ ˜ φ ( t ) , (113)47here we have assumed that ˜ φ ( t ) is already normalized by the energy. We arethus left determining the phase δ such that (108) is satisfied. This is obtainedfor tan δ = − ˜ A i ˜ A r , (114)where ˜ A r and ˜ A i correspond to A r and A i with only φ ( t ) replaced by ˜ φ ( t ). A.2 Infinite domains
As mentioned in section (3), when considering flows in infinite or semi-infiniteintervals, we are not in possession of any theoretical result giving an estimateon the number of discrete eigenvalues. Therefore, we need to carefully check thenumerical results when solving equation (37) by varying the number of basispolynomials and the extend h at which we truncate the numerical domain. Inthe following, we consider VKD modes for the boundary layer flow under asolitary wave at two points in time, 2 t/ Re = − t/ Re = 2. From figure24, we infer that the boundary layer itself has a thickness of around 3 and6 respectively. In figure 25, the magnitude of the first three VKD modes isplotted with respect to the wall normal distance for two different values of thenumerical cut-off parameter h . The solutions are well converged concerning thenumber of polynomials used, ie. 129.The spatial extend of the zeroth VKD-mode seems to be well resolved for both values of h , as it falls off exponentiallyto zero in z . However, the first and second VKD mode do not display such adrop off but fill all the domain given to them, also when doubling h once more(figure not shown). This suggests that for Re = 316, α = 0 .
369 and 2 t/ Re = 2,equation (37) only possesses a single discrete eigenvalue and eigenfunction. Allother eigenvalues and eigenfunctions obtained have to be considered numericalartifacts. This can also be observed for the eigenvalues, cf. figure 26, where only λ is well converged with respect to h , whereas the others gather closer togetherfor increasing h . However, when increasing α , figure 27, we observe that λ is bifurcating from the bundled eigenvalues and becoming positive, indicatingthe existence of a second discrete mode for Re = 316 and 2 t/ Re = 1. As amatter of fact φ and φ in figure 24 resemble more functions typical for thecontinuous spectrum, cf. reference [9], displaying oscillations in the free streamand attenuating inside the boundary layer. For the flow in question, the zerothVKD-mode, φ , obtains its characteristic shape at around 2 t/ Re = − .
3, cf.figure 28. As can be observed from figure 28, for 2 t/ Re = −
1, the eigenfunction φ displays a single antinode over the entire domain similar to φ at 2 t/ Re = 2in figure 24. Gradually, this antinode morphs into the exponential hump withan extend comparable to the boundary layer thickness. Similarly, all othereigenfunctions lose one of their antinode at this point in time (not shown).48 z (a) 2 t/ Re = − −0.20 −0.15 −0.10 −0.05 0.00 0.05 0.10U051015202530 z (b) 2 t/ Re = 2 Figure 24: Velocity profile of the boundary layer under a solitary wave at2 t/ Re = − t/ Re = 2 z |φ ||φ ||φ | (a) h = 30 z |φ ||φ ||φ | (b) h = 120 Figure 25: Magnitude of the first three DVK modes with α = 0 .
369 and β = 0for the boundary layer flow under a solitary wave for Re = 316 at 2 t/ Re = 2.The numerical domain is truncated at different values of h .49 .000.010.020.030.040.050.06−4 −2 0 2 4 6 82t/Re−0.0010−0.0008−0.0006−0.0004−0.00020.0000 m λ λ λ λ (a) h = 30 m λ λ λ λ (b) h = 120 Figure 26: Temporal evolution of characteristic quantities of the nonmodalTollmien-Schlichting wave with α = 0 . β = 0, for the boundary layer flowunder a solitary wave at Re = 316. The numerical domain is truncated atdifferent values of h . λ λ λ λ m (a) h = 30 λ λ λ λ m (b) h = 120 Figure 27: Growth rates λ i and dispersion measure m , max for the boundarylayer flow under a solitary wave in function of α at time 2 t/ Re = 1. Thenumerical domain is truncated at different values of h .50 .00 0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.40 0.45|φ |020406080100120 z Figure 28: Magnitude of the zeroth DVK mode for the boundary layer flowunder a solitary wave for Re = 316 for different values of time.51 Path integral expansion
As touched upon in section 3, equation (69) modeling the nonmodal Tollmien-Schlichting waves as a linear combination of VKD modes, can be considered asystem composed of a base Hamiltonian and a perturbation potential. In thefollowing, we shall assume for simplicity that the base flow is steady and thusmatrices Λ and N are constant in time. Matrix N can be diagonalized: N = SHS T , (115)where H is a diagonal matrix with real eigenvalues on its diagonal. This allowsus to write equation (69) as ˙d = i Hd + Vd , (116)where we have written: d = S T c (117) V = S T ΛS (118)Likewise we have for the fundamental solution X : ˙X = i HX + VX , (119)An integrating factor for equation (119) is given by X = e i H t Y , (120)allowing us to write equation (119) as ˙Y = e − i H t V e i H t | {z } = U ( t ) Y , (121)Path integral formulations, cf. reference [12], allow us to find approximationsfor Y . Integrating equation (119) from t to t , we obtain: Y ( t, t ) = I + t Z t dt U ( t ) Y ( t , t ) (122)Substituting Y back into (121) and performing repeated integrations leads toan expansion formula for Y : Y = I + t Z t dt U ( t ) | {z } A + t Z t dt t Z t dt U ( t ) U ( t ) | {z } B + . . . (123)= ∞ X n =0 t Z t dt . . . t n − Z t dt n U ( t ) . . . U ( t n ) . (124)52he elements of A and B are given explicitly as: A nm = t Z t dt U ( t ) nm (125)= (cid:26) V nn ( t − t ) n = m h m − h n ) V nm (cid:2) e i( h m − h n ) t − e i( h m − h n ) t (cid:3) n = m (126) B nm = t Z t dt t Z t dt U ( t ) U ( t ) nm (127)= P k = n | V nk | ( h n − h k ) (cid:2) h k − h n ) ( t − t ) − e i( h k − h n )( t − t ) (cid:3) + | V nn | ( t − t ) n = m P k = n,m V nk V km ( h m − h k )( h m − h n )( h k − h n ) n ( h k − h m ) e i( h m − h n ) t +( h m − h n ) e i( h k − h n ) t +i( h m − h k ) t + ( h n − h k ) e i( h m − h n ) t o + V nn V nm ( h m − h n ) n i( h m − h n )( t − t ) e i( h m − h n ) t − e i( h m − h n ) t + e i( h m − h n ) t o + V nm V mm ( h m − h n ) n − i( h m − h n )( t − t ) e i( h m − h n ) t + e i( h m − h n ) t − e i( h m − h n ) t o n = m (128)In quantum electrodynamics, the terms in formula (123) represent differentlevels of interaction between the modes of the base Hamiltonian H , equation(115). The first term on the right hand side of (123) represents the undisturbedsolution to the base Hamiltonian. The matrix A , equation (127), on the otherhand, stands for the effect of pairwise interactions between modes under actionof the potential V . Likewise, the matrix B , equation (128), accounts for tripleinteractions between modes under action of the potential V , and so on. Theenergy c † c of the nonmodal Tollmien-Schlichting wave is given by: c † c = d † S † Sd (129)= d † d (130)= d † X † Xd (131)= d † Y † Yd (132)= d † (cid:0) I + A † + B † + . . . (cid:1) ( I + A + B + . . . ) d (133)= d † d + d † (cid:0) A † + A (cid:1) d + d † (cid:0) A † A (cid:1) d + d † (cid:0) B † + B (cid:1) d + d † (cid:0) A † B + B † A (cid:1) d + . . . (134)53 t G GG Figure 29: Amplification G of the optimal perturbation and approximation G ,equation (136) for α = 1 .
21 and β = 0 for Couette flow at Re = 1000.By means of this formula, we can define successive approximations to G : G = 1 (135) G = max d d † (cid:0) I + A † + A (cid:1) d d † d (136) G = max d d † (cid:0) I + A † + A + A † A + B † + B (cid:1) d d † d (137) . . . (138)In figure 29, the amplification G of the optimal perturbation is displayed nextto its approximation G for α = 1 .
21 and β = 0 for Couette flow at Re = 1000.For small times the approximation G follows the solution G closely, but thendiverges for larger times. A typical issue for path integral approximations is thatfor higher orders the approximation diverges. This is also the case for nonmodalTollmien-Schlichting waves, where G diverges even for smaller times than G due to the term A † A . In fact, because of the viscous term, the matrix A haslarge negative eigenvalues which leads to A † A having large positive eigenvaluescausing unphysically large growth. 54 eferences [1] F. Bertolotti, T. Herbert, and P. Spalart. Linear and nonlinear stabilityof the Blasius boundary layer. Journal of Fluid Mechanics , 242:441–474,1992.[2] L. Brandt, P. Schlatter, and D. S. Henningson. Transition in boundary lay-ers subject to free-stream turbulence.
Journal of Fluid Mechanics , 517:167–198, 2004.[3] K. M. Butler and B. F. Farrell. Three-dimensional optimal perturbationsin viscous shear flow.
Physics of Fluids A , 4:1637–1650, 1992.[4] F. Charru.
Hydrodynamic instabilities . Cambridge University Press, 2011.[5] S. D. Chatterji.
Cours d’Analyse 3: Equations diff´erentielles ordinaires etaux d´eriv´ees partielles . Presses polytechniques et universitaires romandes,1998.[6] S. H. Davis and C. von Kerczek. A reformulation of energy stability theory.
Archive for Rational Mechanics and Analysis , 52(2):112–117, June 1973.[7] P. G. Drazin and W. H. Reid.
Hydrodynamic Stability . Cambridge Univer-sity Press, 1981.[8] H. Gustavsson. Energy growth of three-dimensional disturbances in planePoiseuille flow.
Journal of Fluid Mechanics , 224:241–260, 1991.[9] M. J. P. Hack and T. A. Zaki. The continuous spectrum of time-harmonicshear layers.
Physics of Fluids , 24:034101–1–23, 2012.[10] J. Jimenez. How linear is wall-bounded turbulence?
Physics of Fluids ,25:110814–1–19, 2013.[11] P. Luchini and A. Bottaro. Adjoint equations in stability analysis.
AnnualReview of Fluid Mechanics , 46:493–517, 2014.[12] P. A. Martin and F. Rothen.
Probl`emes `a N-corps et champs quantiques .Presses polytechniques et universitaires romandes, 1990.[13] P. J. Schmid. Nonmodal stability theory.
Annual Review of Fluid Mechan-ics , 39:129–162, 2007.[14] P. J. Schmid and D. S. Henningson.
Stability and Transition in Shear Flows .Springer Science+Business Media, 2001.[15] J. Shen. Efficient spectral-galerkin method i. direct solvers for the secondand fourth order equations using legendre polynomials.
Siam Journal ofScientific Coputing , 15:1489–1505, 1994.5516] J. Shen. Efficient spectral-galerkin method ii. direct solvers of second fourthorder equations by using chebyshev polynomials.
SIAM Journal of Scien-tific Computing , 16(1):74–87, 1995.[17] B. M. Sumer, P. M. Jensen, L. B. Sørensen, J. Fredsøe, P. L.-F. Liu, andS. Carstensen. Coherent structures in wave boundary layers. part 2. solitarymotion.
Journal of Fluid Mechanics , 646:207–231, 2010.[18] L. N. Trefethen, A. E. Trefethen, S. C. Reddy, and T. A. Driscoll. Hydro-dynamic stability without eigenvalues.
Science , 261:578–584, 1993.[19] E. Ventsel and T. Krauthammer.
Thin Plates and Shells . Marcel Dekker,2001.[20] J. C. G. Verschaeve, G. K. Pedersen, and C. Tropea. Non-modal stabil-ity analysis of the boundary layer under solitary waves.