Boundary-layer instability is influenced by non-parallel-flow effects, that is, the streamwise variation of the mean flow not only modifies directly the growth rate but also distorts the shape (i.e., the transverse distribution) of the disturbance, thereby affecting the growth rate indirectly. In this paper, we present a simple but effective local approach to spatial instability of three-dimensional boundary layers, which takes into account both direct and indirect effects of non-parallelism. Unlike existing WKB/multi-scale type of methods, the present approach is non-perturbative in which non-parallelism does not need to be a higher-order correction to the leading-order prediction by the parallel-flow approximation. The non-parallel-flow effects are accounted for by expanding the local mean flow and the shape function as Taylor series, and this leads to a sequence of extended eigenvalue problems, depending on the order of truncation of the Taylor series. These eigenvalue problems can be solved effectively by standard numerical methods. In the case of the Blasius boundary layer, the predictions are verified and confirmed by direct numerical simulation results as well as by the non-parallel theory of Gaster. The non-parallel-flow effects on the eigenvalues, eigenfunctions, and neutral curves for planar and oblique Tollmien-Schlichting (T-S) waves are discussed. The distortion of the eigenfunction is found to have a significant effect, which may be stabilizing or destabilizing depending on the ranges of the Reynolds number and frequency.

The problem of instability and laminar-turbulent transition of boundary layers has been studied extensively for many decades because of its fundamental importance in classical physics and practical relevance to many modern technologies, such as designs of gas-turbine-engine blades, vanes, and aircraft wings. While much progress has been made towards understanding its many aspects, the problem as a whole has not been fully resolved and remains a topic of continuing interest.

The problem proves challenging because it involves several processes occurring on different length and time scales. When a uniform flow passes an obstacle such as an airfoil, a thin shear (boundary) layer forms near the surface, and its thickness increases with the downstream distance x. At a distance x, the thickness is of order Re−1/2x, much smaller than x when the Reynolds number Re = U x / ν is large, where U is the free-stream velocity and ν the kinematic viscosity. It is understood that transition is controlled by both external disturbances (such as free-stream turbulence) and the inherent instability of the mean flow, the boundary layer. In a low-turbulence environment, transition starts with external disturbances exciting instability waves, whose length scale is much shorter than x (but much longer than the boundary-layer thickness). These instability waves amplify exponentially downstream, then undergo nonlinear evolution, and finally break down, through rapid secondary instability, into small-scale random motions.

Instability of the boundary-layer is an important aspect of transition. Since disturbances at early stage are pretty weak, linear stability theory (LST) has been developed.1 The basic idea behind LST is to introduce into the mean flow a small perturbation. Substitution into the Naiver-Stokes (N-S) equations followed by linearization about the base state leads to the linearized N-S equations governing the disturbance. Stability/instability is then delineated according to whether perturbations decay or grow. In general, this requires solving an initial-value problem, which is difficult and inconvenient. For parallel flows, one may seek so-called normal-mode solutions, or equivalently take Fourier transform with respect to the streamwise and spanwise variables, thereby reducing the linearized N-S equations into a set of ordinary differential equations, which are further simplified to the well-known Orr-Sommerfeld (O-S) equation. The latter forms, along with the homogeneous boundary conditions, an eigenvalue problem. Boundary layers are unfortunately non-parallel flows. The eigenvalue problem is tenable only by making the local parallel-flow assumption, under which the streamwise dependence of the mean flow is taken to be parametric, and the streamwise variation and the normal velocity are ignored. Typically, spatial instability is formulated, in which the perturbation is assumed to be time periodic but evolves downstream. This leads to an O-S equation for a streamwise velocity profile u ̄ frozen at each location. The eigenvalue problem is solved to determine the local growth/decay rates of eigen modes, referred to as Tollmien-Schlichting (T-S) waves in the case of the flat-plate boundary layer.2 Integrating the local growth rate with respect to the streamwise variable yields a measure of the accumulated amplification, or the so-called N-factor, which has been used to estimate transition location on an empirical basis.

The parallel-flow approximation usually gives qualitatively correct predictions. For the incompressible flat-plate boundary layer, the predicted neutral curve for planar modes agrees well quantitatively with the experiment data for Re1/2 > 600,3,4 but appreciable discrepancy exists for Re1/2 < 600 due to non-parallel-flow effects.1 It is recognized that non-parallelism is more significant for oblique T-S waves5,6 and cross-flow vortices.7 

Non-parallel-flow effects on the growth rate consist of two parts: (a) the direct contribution of the streamwise variation and the normal velocity of the mean flow, u ̄ / x and v ̄ and (b) the influence of the streamwise distortion of the local eigenfunction. The contribution (a) can easily be accounted for by including u ̄ / x and v ̄ in the perturbation equation, leading to the extended O-S equation (EOS).8 The solution suggested a destabilizing role in the sense that the lower branch of the neutral curve shifts to lower Reynolds numbers and the unstable region becomes broader. Several different approaches have been developed to take both effects (a) and (b) into account. The first approach is based on the technique of WKB/multi-scale expansions; see, for example, Boutheir,9 Ling and Reynolds,10 Gaster,11 Saric and Nayfeh,12 Stijn and Vooren,13 and Thomas and Philip.14 While these analyses were presented and/or interpreted somewhat differently, they were all based on the common idea that the background mean-flow is varying slowly on a length scale much longer than the characteristic wavelength of the instability modes, (which were implicitly assumed to be comparable with the boundary-layer thickness, a point to which we shall return). Assume that the ratio of the length scales is ϵ−1, with ϵ ≪ 1. It follows that mathematically a WKB type of approximation is applicable, according to which the instability is governed, at leading order, by the usual O-S equation in the parallel-flow theory. The non-parallel correction is of O(ϵ) and is governed by an inhomogeneous O-S equation; the solvability condition of which allows for determination of either the correction to the local growth rate or the shift of the neutral position for a mode with a given frequency. The latter was calculated by Ling and Reynolds,10 who applied this idea specifically for the parameters (i.e., the frequency and Reynolds number) in the vicinity of the neutral curve. In addition, by expanding the mean flow as a Taylor series about the leading-order approximation of the neutral point, they also obtained the local variation of the growth rate. Gaster11 generalized the multi-scale analysis at an arbitrary streamwise location (rather than being restricted to the neighbourhood of the neutral curve) by introducing an amplitude function of the slow variable ξϵx to account for the non-parallel-flow effects. Gaster pointed out that the neutral curve depends on which the quantity is used to calculate the local growth rate. For large Reynolds numbers, the parallel-flow approximation was found to be accurate enough, but the prediction by the non-parallel theory agrees better with experiments at moderate Reynolds numbers. Later, direct numerical simulations (DNSs)15 and calculations using parabolized stability equations (PSE)7 confirmed that the theory of Gaster11 was remarkably accurate (see below). Saric and Nayfeh12 and Crighton and Gaster16 applied the WKB/multi-scale approach to boundary layers with pressure gradients and to free shear layers.

In all the analyses above, the variation of the mean flow is assumed to be slow on the basis that ϵ ≪ 1, while at the same time, the seemingly O(Re−1/2) viscous effect is retained at leading order, that is, ϵ and Re−1/2 are artificially treated as being independent parameters despite that ϵ = O(Re−1/2). If for asymptotic consistency the viscous terms are neglected as well, the resulting inviscid solution would become invalid near the wall and also in the so-called critical layer, where viscosity plays a leading-order role. A multi-deck vertical structure then arises, which may take different forms depending on whether lower or upper branch modes are considered. For lower-branch modes, the characteristic wavelength is of order Re1/8 times the boundary-layer thickness, and the well-known triple-deck structure emerges.17,18 Using this formalism, Smith18 constructed the high-order asymptotic solution to the linearized N-S equations. His analysis shows that the correction by non-parallelism to the growth rate is of O(Re−3/8), apparently much larger than O(Re−1/2) that the method of Bouthier9 and Gaster11 suggests. This seemingly inconsistent result arises because the characteristic wavelength is taken to be comparable with the boundary-layer thickness in the latter method, but it does not mean that the prediction by the latter approach was wrong since once a correct scaling of the wavelength is used in the result, the correct and consistent order of magnitude for the non-parallelism follows. Upper-branch modes are governed by an asymptotic structure consisting of five decks. By systematically constructing high-order expansions in each of these zones and performing matching, Bodonyi and Smith19 calculated the non-parallel-flow correction.

While the multi-layered asymptotic formalism is consistent and systematic, the analysis, involving asymptotic matching between different decks in the normal direction, is extremely complex, and the prediction, in either of the lower- and upper-branch regimes, is not accurate at moderate Reynolds numbers of practical relevance. In contrast, the WKB/multi-scale method retains tacitly the viscous effect in the entire boundary layer so that the resulting O-S equation serves as a uniformly valid composite approximation, alleviating the need of analyzing different vertical zones. As far as linear instability is concerned, this method has been shown to predict non-parallelism rather satisfactorily. However, it should be pointed out that if nonlinear effects are to be considered, the high-Reynolds-number asymptotic approach is preferred because the WKB/multi-scale method masks the precise nature of the eigenfunction in the critical layer and viscous wall layer, and as a consequence fails to account for the delicate balance between inertia, viscosity and nonlinearity.

Bertolotti et al.20 proposed a local approach that accounts for non-parallelism. The mean flow and the solution for the disturbance are expanded as Taylor series of the variable xx0, for example, the disturbance shape function ϕ(x, y) = ϕ0(y) + (xx0) ϕ1(y) + O((xx0)2), as was done in Ling and Reynolds,10 but the main differences are that x0 is arbitrary rather being a first approximation of the neutral position and no asymptotic expansion with respect to the small parameter ϵ ≪ 1 was performed at outset. By collecting and balancing terms of order O(1) and O[(xx0)], they obtained two coupled equations for ϕ0 and ϕ1, which form an eigenvalue problem with (ϕ0, ϕ1) appearing as the eigenfunction. This is then solved by an iterative procedure based on the fact that the parameter ϵ ≪ 1. The end result is essentially equivalent to that given by the WKB/multi-scale method.

The approach of parabolized stability equations (PSE)7 has become a popular tool for analyzing instability of weakly non-parallel flows such as boundary layers. In this approach, the fast streamwise variation of a disturbance is assumed to be sinusoidal and characterised by a local wavenumber so that its solution is written as a product of a carrier wave and a shape function varying slowly in the streamwise direction. This is similar to the WKB/multi-scale method, but the difference is that no explicit reference is made to the small parameter ϵ measuring the ratio of the length scales, and the solution is not sought via an asymptotic expansion. Rather, substituting the assumed form into the N-S equations and neglecting the second-order derivative of the shape function with respect to the streamwise variable, one obtains the perturbation equations, which are parabolic and hence can be solved effectively by a marching method. But these equations are not closed because the wavenumber is not yet known. A heuristic normalization operation is used to define a local wavenumber at each streamwise position. In order to start the marching procedure, an initial condition must be specified. It is sometimes claimed that the method of PSE is applicable to any form of initial conditions. However, if an arbitrary initial disturbance is imposed, the evolution experiences a rapid transient phase before approaching a history-independent asymptotic state. The history-dependent transient behaviour is most likely to be meaningless because the rapid adjustment violates the assumption that the shape is slowly varying. In most applications, the initial condition is chosen to be local eigen modes. However, if their wavenumbers and shape functions are given by LST with non-parallel effects being ignored, the disturbance would also exhibit rapid adjustment near the inlet unless the Reynolds number is extremely large. In order to ensure a gradual development, consistent with the required assumption of the method, it is necessary to prescribe an inlet condition which accounts for effects of non-parallelism. This can be done using the WKB/multi-scale approach or the local procedure of Bertolotti et al.20 In this case, PSE reproduce what local non-parallel theories predict. However, one may argue that since a history-independent state appears to exist, a method capable of computing it directly is preferred. Turner and Hammerton21 demonstrate that for relatively low-frequency modes, the upstream Lam-Rott eigen solutions, which account for non-parallelism, match with the PSE solution in an overlapping region and may therefore serve as the appropriate initial condition.

DNSs, which solve the complete N-S equations, have been performed by Fasel and Konzelmann15 to investigate non-parallel-flow effects. They found that the growth rate and the neutral curve depend strongly on the chosen criteria, and that the deviation of DNS and theoretical results from experiment measurements cannot be all attributed to non-parallelism; other factors, such as small pressure gradient, may be responsible too. Among available non-parallel instability theories, that of Gaster11 agrees remarkably well with the DNS results. While DNS stands as a powerful tool for studying instability of boundary layers, the computation time is long and the cost is high and so, it is usually used to assess the validity of theoretical results rather than for comprehensive parametric studies.

In this paper, we present a simple approach for studying the non-parallel-flow effects on spatial instability of shear flows. Unlike the WKB/multi-scale type of method, the approach is non-perturbative in which non-parallelism does not need to be a small correction to the leading-order prediction by the parallel-flow approximation. Such situations do exist. A recent analysis by Butler and Wu,30 based on asymptotically large Reynolds numbers, shows that the growth rates of lower-branch stationary cross-flow vortices in three-dimensional boundary layers are influenced at leading order by non-parallelism. Furthermore, the non-parallelism is found to be associated with the distortion of the eigenfunction. That effect must persist and may be substantial at finite Reynolds numbers but cannot possibly be accounted for by an extended O-S equation. The present theory can be applied to investigate such leading-order effects of non-parallelism. Compared with the PSE approach, the method is local in which it computes directly the growth rates of asymptotic modes.

In our method, all terms associated with the streamwise variation of the mean flow are taken into account, whereas only certain such terms are included in existing approaches. The shape function of the disturbance, ϕ, is allowed to be dependent on both the normal and streamwise coordinates, x and y. Its local variation in the vicinity of an arbitrary point xa with x is then approximated by a Taylor expansion,

ϕ ( x , y ) = ϕ 0 ( y ) + ( x x a ) ϕ 1 ( y ) + ( x x a ) 2 ϕ 2 ( y ) + .

Substituting the solution of the assumed form into the linearized N-S equations leads to an extended eigenvalue (EEV) problem with the extended vector (ϕ0, ϕ1, ϕ2) appearing as the eigenvector. This is solved directly, and the effective local wavenumber and growth rate, with the effects of non-parallelism being accounted for, are calculated immediately. The method will therefore be referred to as “EEV,” approach. It is simpler than the perturbative WKB/multi-scale method because it avoids calculating the adjoint function and evaluating various integrals.

The formulation will be presented for a general three-dimensional boundary layer in Sec. II. We show that accounting for the shape distortion of the disturbance by a local expansion leads to an extended eigenvalue problem. In Sec. III, numerical issues concerning the eigenvalue problems and DNS, which we will perform to verify the theoretical results, are highlighted. Calculations will be carried out for the Blasius boundary layer, for which experimental and DNS data are available. In Sec. IV, numerical solutions will be presented, and non-parallel-flow effects on the eigenvalues, eigenfunctions, and neutral curve are investigated for planar and oblique T-S waves. Comparisons will be made with the prediction by the non-parallel theory of Gaster,11 with the DNS results of Fasel and Konzelmann15 as well as with the experimental data of Kachanov et al.5,6 Concluding remarks are given in Sec. V, where we indicate further applications of the present method.

The concern of this paper is with spatial instability of weakly non-parallel shear flows. To fix the idea, we consider the three-dimensional incompressible boundary layer that develops over a surface. The flow is described in orthogonal curvilinear coordinates (x, y, z), in which the surface is located at y = 0 with its leading edge at x = 0, where the x axis points to the downstream direction. The spanwise coordinate is z. Let U , ρ , and μ denote the free-scream velocity, density, and viscous coefficient, respectively, where the superscript ∗ indicates a dimensional quantity. The spatial and time variables will be non-dimensionalized by δ 0 and δ 0 / U , respectively, where δ 0 is the displacement thickness of the boundary layer at a typical location x 0 ; the resulting dimensionless coordinates and time will be denoted as (x, y, z) and t, respectively. The velocity and pressure are normalized by U and ρ U 2 respectively, and the corresponding non-dimensional quantities are denoted by (u, v, w) and p.

The flow is characterized by a Reynolds number R, defined as

R = ρ U δ 0 / μ .
(1)

The dimensionless three-dimensional incompressible N-S equations read

u x + v y + w z = 0 ,
(2)
u t + u u x + v u y + w u z + p x = 1 R ( 2 u x 2 + 2 u y 2 + 2 u z 2 ) ,
(3)
v t + u v x + v v y + w v z + p y = 1 R ( 2 v x 2 + 2 v y 2 + 2 v z 2 ) ,
(4)
w t + u w x + v w y + w w z + p z = 1 R ( 2 w x 2 + 2 w y 2 + 2 w z 2 ) .
(5)

For simplicity, here, we restrict to boundary layers over a flat surface, but the formulation and method to be presented can easily be generalized to include curvature effects by using the orthogonal curvilinear coordinates with appropriate Lamé coefficients.

Let the three-dimensional mean flow field be denoted by Q ̄ ( x , y ) = ( u ̄ , v ̄ , w ̄ , p ̄ ) . In order to study its stability, the mean flow is perturbed by a three-dimensional small-amplitude disturbance Q′(x, y, z, t). The total flow field Q(x, y, z, t) is written as

( u , v , w , p ) = ( u ̄ ( x , y ) , v ̄ ( x , y ) , w ̄ ( x , y ) , p ̄ ( x , y ) ) + ( u , v , w , p ) Q ̄ ( x , y ) + Q ( x , y , z , t ) .
(6)

Substituting Eq. (6) into Eqs. (2)(5) and linearizing about the mean flow, we obtain the linear N-S equations for the perturbation. For ensuing analysis, it is convenient to cast the perturbation equations into the matrix form

[ A + B 1 x + B 2 y + B 3 z + C ( R t + 2 x 2 + 2 y 2 + 2 z 2 ) ] Q = 0 ,
(7)

here we have put

A ( x , y ) = 0 0 0 0 u ̄ / x u ̄ / y 0 0 v ̄ / x v ̄ / y 0 0 w ̄ / x w ̄ / y 0 0 , B 1 ( x , y ) = 1 0 0 0 u ̄ 0 0 1 0 u ̄ 0 0 0 0 u ̄ 0 ,
(8)
B 2 ( x , y ) = 0 1 0 0 v ̄ 0 0 0 0 v ̄ 0 1 0 0 v ̄ 0 , B 3 ( x , y ) = 0 0 1 0 w ̄ 0 0 0 0 w ̄ 0 0 0 0 w ̄ 1 , C = 1 R 0 0 0 0 1 0 0 0 0 1 0 0 0 0 1 0 .
(9)

In LST for parallel flows, the disturbance quantity Q′(x, y, z, t) is assumed to take the normal mode form ϕ(y)exp{i(αx + βzωt)} with the shape function ϕ depending on y only, where the frequency, ω, and the streamwise and spanwise wavenumbers, α and β, are all constants. For a non-parallel flow, the mean flow varies with x (associated with which a vertical velocity component v ̄ is also present). The perturbation Q′(x, y, z, t) can be written as

Q ( x , y , z , t ) = ϕ ( x , y ) exp { i ( x α ( x ) d x + β z ω t ) } + c.c. ,
(10)

where c.c. stands for complex conjugate, ω and β are constants, but the streamwise wavenumber α and the shape function ϕ(x, y) depend on x.

Substitution of Eq. (10) into Eq. (7) yields the linear stability equations

( A ˆ p + A ˆ n ) ϕ + B ˆ p ϕ x + C ˆ p ( 2 ϕ x 2 + α x ϕ ) = 0 ,
(11)

where

(12)
ϕ ( x , y ) = u ˆ v ˆ w ˆ p ˆ , A ˆ p ( x , y ) = i α D i β 0 s u ̄ y 0 i α 0 s 0 D 0 w ̄ y s i β ,
A ˆ n ( x , y ) = 0 0 0 0 u ̄ x + v ̄ D 0 0 0 v ̄ x v ̄ y + v ̄ D 0 0 w ̄ x 0 v ̄ D 0 ,

B ˆ p ( x , y ) = 1 0 0 0 u ̄ 2 i α / R 0 0 1 0 u ̄ 2 i α / R 0 0 0 0 u ̄ 2 i α / R 0 , C ˆ p = 1 R 0 0 0 0 1 0 0 0 0 1 0 0 0 0 1 0 ,
(13)

here we have put

s = i α u ̄ + i β w ̄ i ω [ ( i α ) 2 + D 2 + ( i β ) 2 ] / R , D = / y .

The matrices A ˆ p , B ˆ p , and C ˆ p represent the contribution of the parallel part of the mean flow. The non-parallel effects consist of two parts: the direct contribution of the streamwise variation of the mean flow, represented by the matrix A ˆ n , and the streamwise distortion of the shape eigenfunction, represented by the derivatives of ϕ with respective to x. The streamwise variation of α can be absorbed into the shape function as will be shown later.

For a parallel mean flow, or under the local parallel-flow assumption, the streamwise variations of the mean flow, the shape function ϕ, and the wavenumber α are absent or ignored, and consequently, the linear stability equation (11) simplifies to

[ A ˆ p ] x a ϕ ( x a , y ) = 0 ,
(14)

which holds at an arbitrary location xa and can be further reduced to an ordinary differential equation, namely, the classical O-S equation. Equation (14) with the usual homogeneous boundary conditions poses an eigenvalue problem. The eigenvalue problem is to find the complex eigenvalue α(xa) and the corresponding eigenvector ϕ(xa, y) at a location xa for a given β and ω. The eigenvector ϕ(xa, y) is usually normalized by u ˆ ( x a , y m ) , where ym is the location at which | u ˆ | reaches its maximum. This is what linear stability analysis is about, and the prediction will be referred to (or labelled) as LST in this paper.

In the case of a non-parallel flow, the streamwise variations of both the mean flow and the shape function are taken into account in the linear stability Equation (11). Unfortunately, this set of partial differential equations is elliptic and is extremely difficulty to solve. In order to simplify the problem, the local variation with respect to x in the vicinity of an arbitrary point xa is approximated by a Taylor expansion, namely,

ϕ ( x , y ) exp [ i x α ( x ) d x ] = { ϕ 0 ( y ) + ( x x a ) ϕ 1 ( y ) + ( x x a ) 2 2 ϕ 2 ( y ) + } exp [ i x α ( x a ) d x ] ,
(15)

in which the streamwise variation of α has been absorbed into the O((xxa)2) term in the expansion of the shape function. It follows that ∂α/∂x in Eq. (11) can be set to zero since α = α(xa) is a constant in the vicinity of xa. This term is also of smaller size than other terms representing non-parallelism when the Reynolds number is large.

Similarly, the variation of the mean flow is accounted for by expanding its velocities, u ̄ , v ̄ , and w ̄ , as Taylor series. Correspondingly, the operator expands as

A ˆ p ( x ) = A ˆ p ( x a ) + ( x x a ) A ˆ p x + ( x x a ) 2 2 2 A ˆ p x 2 + ,
(16)

and similar expansions are performed for B ˆ p and A ˆ n . Introducing the approximations of (15) and (16) into Eq. (11) and retaining the first a few terms, we obtain a sequence of extended eigenvalue problems. The method and the solution will be referred to as EEVn, where n refers to the order at which the Taylor series are truncated. The idea of the local expansion has its root in the work of Ling and Reynolds10 and Bertolotti et al.,20 but it is now being developed systematically, and extended to three-dimensional disturbances and three-dimensional mean flows. The extension to higher orders would allow for predictions with increasing accuracy.

1. EOS (EEV0)

If the streamwise variation of the mean flow is included, but the shape distortion is discarded, i.e., only the first term on the right-hand side of Eq. (15) is retained, then the eigenvalue problem reads

[ A ˆ p + A ˆ n ] x a ϕ 0 = 0 .
(17)

This kind of approximation was first introduced by Barry and Ross.8 As the system Eq. (17) can be reduced to an extended O-S equation, the method and solution will be referred to as EOS, but may alternatively be designated as EEV0 in the present framework. It should be pointed out that EOS (or EEV0) is not a complete or consistent approximation because the neglected effect of shape distortion has, from asymptotic viewpoint, the same order of magnitude as the direct influence, and numerically turns to be greater, as will be shown later. The prediction of EOS will be presented only for the purpose of illustrating the importance of shape distortion.

2. EEV1

When the non-parallel contribution of the mean flow and the first two terms on the right-hand side of Eq. (15) are retained, we obtain an extended eigenvalue problem,

A ˆ p + A ˆ n B ˆ p ( A ˆ p + A ˆ n ) x A ˆ p + A ˆ n + B ˆ p x x a ϕ 0 ϕ 1 = 0 0 ,
(18)

in which the extended vector (ϕ0, ϕ1)T appears as the eigenfunction. At the present level of approximation, the shape distortion corresponding to ∂ϕ/∂x is accounted for, but ∂2ϕ/∂x2 in Eq. (11) is set to zero. Similar to LST, the extended eigenvalue problem can be solved directly with the mean flow information at the location x = xa.

Bertolotti et al.20 neglected some of the terms representing the non-parallel mean flow and obtained a similar equation

A ˆ p + A ˆ n B ˆ p A ˆ p x A ˆ p x a ϕ 0 ϕ 1 = 0 0 .
(19)

In their work, an iterative procedure was employed to solve the eigenvalue problem Eq. (19). In the present work, all terms representing non-parallelism are retained in the spirit of developing a non-perturbative approach, and the eigenvalue problem Eq. (18) will be solved directly. For comparison purpose, we will also solve Eq. (19).

3. EEV2

We now retain the non-parallel contribution of the mean flow and the first three terms on the right-hand side of Eq. (15). There follows an eigenvalue problem

A ˆ p + A ˆ n B ˆ p C ˆ p ( A ˆ p + A ˆ n ) x A ˆ p + A ˆ n + B ˆ p x B ˆ p 2 ( A ˆ p + A ˆ n ) x 2 2 ( A ˆ p + A ˆ n ) x + 2 B ˆ p x 2 A ˆ p + A ˆ n + 2 B ˆ p x x a ϕ 0 ϕ 1 ϕ 2 = 0 0 0 ,
(20)

in which the extended vector (ϕ0, ϕ1, ϕ2)T appears as the eigenfunction. The approximation at this level accounts for the shape distortion corresponding to ∂ϕ/∂x and ∂2ϕ/∂x2. Similar to LST and EEV1, this eigenvalue problem can be solved directly with the local information of the mean flow at x = xa.

The homogeneous systems, Eqs. (14), (17), (18), and (20), are all subject to the usual boundary conditions

( u ˆ , v ˆ , w ˆ ) = ( 0 , 0 , 0 ) at y = 0 , ( u ˆ , v ˆ , w ˆ , p ˆ ) ( 0 , 0 , 0 , 0 ) as y ,
(21)

which correspond to no-slip on the wall and vanishing of the perturbation at infinity, respectively. Among the resulting eigenvalue problems, the first will be is the standard LST, while the remaining ones will be referred to as EOS, EVV1, and EVV2, respectively.

In principle, it is possible to extend the formulation to an arbitrary order, leading to methods capable of making prediction with arbitrarily small accuracy. However, EEV1 and EEV2 usually suffice for moderate Reynolds numbers of interest, as will be seen later.

The present formulation is obviously rather general, and is applicable, with minimal adjustment, to general weakly non-parallel three-dimensional shear flows, such as free shear layers, in which case, the only modification is to replace the first boundary condition in (21) by ( u ˆ , v ˆ , w ˆ , p ˆ ) ( 0 , 0 , 0 , 0 ) as y → −∞.

For fixed values of the Reynolds number R, the frequency ω, and the spanwise wavenumber β, the eigenvalue α and its corresponding eigenfunction can be determined by solving one of those eigenvalue problems. By changing the location xa, each eigenvalue problem can be solved for the whole domain.

The growth rate of a physical disturbance quantity q′ is defined as its logarithmic derivative,

1 q q x = ln q x .
(22)

For a parallel flow, all disturbance quantities grow or decay at exactly the same rate, which is −αi, the imaginary part of the eigenvalue of the O-S equation. For a non-parallel mean flow, the growth rate depends on the eigenvalue α, the choice of the quantity q′, and the location y. Since q = q ˆ exp { x α d x } , once α and q′ are found by solving the extended eigenvalue problem, the corrected wavenumber is given by

α ̃ ( y ) = α + α ˆ ( y ) , α ˆ ( y ) = 1 q ˆ ( y ) q ˆ ( y ) i x ,
(23)

where α ˆ stands for the correction to the wavenumber due to the distortion of the shape function q ˆ ( y ) . Usually, the streamwise velocity u ˆ of the perturbation is chosen and y is taken to be the peak position of | u ˆ | since u ˆ is the quantity that can be measured rather accurately. Unless stated otherwise, this is how the corrected wavenumber α ̃ will be calculated.

Note that rather than being treated as a small correction as in WKB/multi-scale approach, in the present method, the effects of non-parallelism, especially that of shape distortion, are accounted for simultaneously along with the leading-order growth rate. Owing to this non-perturbation nature, the method is applicable to problems where non-parallelism plays a leading-order role.

An instability mode is usually referred to by its non-dimensional frequency F, which is related to ω (used in the formulation) by

F = ω / R × 1 0 6 .

In addition to its growth rate α ̃ i , other characteristics of interest are its propagation angle σ and the streamwise phase velocity cx, which are defined, respectively, as

σ = arctan ( β / α ̃ r ) , c x = ω / α ̃ r .
(24)

While the formulation was presented for a three-dimensional boundary layer, the calculations will be performed for the Blasius boundary layer, for which there exists a similarity solution.

The mean flow Q ̄ ( x , y ) , its derivatives with respect to x and y, and hence the matrices appearing in the eigenvalue problems, Eqs. (14), (17), (18), and (20), can be evaluated at an arbitrary location x = xa either directly or by using a five-point finite-difference scheme. A similar scheme is used to discretize the operator D = ∂/∂y in the matrices.

The eigenvalue problems are discretized using also the five-point finite-difference scheme. The absence of the boundary condition for p ˆ at the wall is, as usual, remedied by discretizing the normal pressure gradient, p ˆ / y , on a staggered grid, i.e., at the centre of each interval between two adjacent mesh points. The eigenvalues to the discretized systems are found by Muller’s method, in which the block Gaussian elimination with partial pivoting was used to calculate the determinant, and the latter was made to vanish through an interactive procedure.

In our calculations, the mesh in the y-direction is non-uniform, being gradually stretched with the distance according to

y j = y max 4 η 4 + ( 1 + 8 k 1 3 ) ( 1 η 2 ) , η = j J k 2 1 π ( k 2 + 1 ) sin ( j π / J ) , ( 0 j J ) ,
(25)

where k1 and k2 are parameters. The default grid parameters in Eq. (25) are selected as k1 = 100, k2 = 60, the position of the outer boundary ymax = 100, and the number of points in the y − direction is J = 201. With these parameters, there are more than 120 points in one boundary-layer thickness and the mesh size Δy at the wall is smaller than 0.0025. Tests of grid dependency show that the relative error in the eigenvalue is less than 0.4%, indicating that the mesh size is fine enough.

Numerical calculations indicate that the extended eigenvalue problem EEVn (n = 1, 2) exhibits multiplicity in which it has (n + 1) number of different eigenvalues α with −αi being positive or slightly negative, (which are potentially of physical relevance). However, they lead to the corrected wavenumbers α ̃ , which are identical to graphical precision in EEV1 (for planar modes) and EEV2 (for oblique modes), consistent with the fact that the physical solution of the eigenvalue problem for the Blasius boundary layer ought to be unique. In practice, it suffices to select the eigenvalue α with the largest −αi. A further discussion of multiplicity and uniqueness is given in the appendix.

In the DNS that we performed to verify the EEV approach, the free-stream parameters are selected under the standard atmospheric pressure near the ground, namely, U = 34 m / s , ρ = 1 . 225 kg / m 3 , p = 101 325 Pa , μ = 1 . 79 × 1 0 5 Pa s , and c = 340 m / s , which give the unit Reynolds number Re = 1.1634 × 106/m and a Mach number M = 0.1. For a given Reynolds number R = Re δ 0 , there is a corresponding displacement thickness δ 0 and the streamwise x 0 , which is the inlet location of the computation domain.

The complete N-S equations were solved numerically using the method in Huang et al.22 The convective terms were split and approximated by a fifth-order weak upwind finite-difference scheme, while the viscous terms were discretized by applying twice an eighth-order central difference scheme for the first-order derivative. The governing equations were integrated fully explicitly in time by using a third-order Total Variation Diminishing (TVD) Runge-Kutta method. As the Blasius similarity solution for the mean flow is not a steady solution of the N-S equations, a small steady source term of O(10−6) was added so that the mean flow is a steady state of the modified equations.23 The evolution of the disturbance is equivalent to that in a true mean flow. A T-S wave with a small amplitude A0 = 10−5 and frequency ω (F) was imposed at the inlet of the computation domain. In the x-direction, the grid is uniform with mesh points being xi = x0 + iΔx (0 ≤ iI), where Δx = (xex0)/I with x0 and xe denoting the positions of the inlet and outlet of the computation domain, respectively, and I being the number of mesh points. The mesh size Δx was chosen to make sure that there are at least 70 points within one wavelength. The mesh distribution in the y-direction is the same as Eq. (25). The time step is Δt = 10−3, giving a Courant-Friedrichs-Lewy (CFL) number Cfl = 0.8. For a planar T-S wave, DNS using 80 threads (2.6 GHz) takes 0.658 s for one-hundred time steps, but lasts nearly 4 h for the wave to propagate from the inlet to the outlet. For a three-dimensional mode, the simulations take much longer (about 10 days). In comparison, solving the LST/EEVn equations to obtain the growth rates in the same streamwise (Reynolds number) range requires only about 1 min on one CPU.

In order to validate our computational code and to verify the validity of the EEV approach, we first carried out calculations for EOS (EEV0) and EEV1 formulations. The respective results for the growth rate are compared in Figure 1(a) with those predicted by LST, DNS,15 and the non-parallel theory.11 The EEV1 prediction is in good agreement with the DNS result of Fasel and Konzelmann15 and the theoretical result of Gaster;11 the latter itself has been confirmed by experiments,3,4 DNS15 and PSE.7 However, there is an appreciable difference between EOS and DNS. Despite that EOS has accounted for the direct effect of the streamwise variation of the mean flow, its prediction appears to be closer to LST than to DNS, suggesting that the shape distortion constitutes the dominant part of non-parallel-flow effects.

FIG. 1.

(a) Comparison of the growth rates α ̃ i obtained by EOS and EEV1 with the non-parallel theory of Gaster11 and the DNS of Fasel and Konzelmann;15 the frequency F = 140. (b) Comparison of the growth rates α ̃ i obtained by EEV1 and EEV2 with our own DNS; the frequency F = 60. The inset in (b) is an enlarged view.

FIG. 1.

(a) Comparison of the growth rates α ̃ i obtained by EOS and EEV1 with the non-parallel theory of Gaster11 and the DNS of Fasel and Konzelmann;15 the frequency F = 140. (b) Comparison of the growth rates α ̃ i obtained by EEV1 and EEV2 with our own DNS; the frequency F = 60. The inset in (b) is an enlarged view.

Close modal

In Figure 1(b), we compare the predictions by LST, EEV1, and EEV2 with our own DNS for the frequency F = 60, which is lower than the most unstable frequency F = 96. LST differs appreciably from DNS, but the results from EEV1 and EEV2 are almost identical to DNS, indicating that the non-parallel-flow effects are captured by EEVn (n = 1, 2). Note that we have traced the unstable mode to the subcritical regime upstream, where the mode is highly damped, and the agreement with DNS is good. Farther upstream, mode crossing occurs, that is, the mode being traced is no longer least stable. The eigenvalue calculation then becomes intricate. At the same time, in DNS, the disturbance consists of the newly emerged least damped mode as well as the mode that is not being traced, and an accurate extract of the growth rate of the latter by DNS becomes an issue due to modal interference.

Table I displays the results obtained by different methods, in which the result of LST by Jordinson24 and the local solution by Bertolotti20 is listed as a reference. When R = 998, the growth rate α ̃ i obtained by EOS is 7% larger compared to that by LST, and the EEV1 result is 12% larger than that of LST. There is an appreciable difference between EEV1 and EOS, indicating that the effect of shape distortion is significant. As the Reynolds number R increases, the difference between the LST and EEV1 predictions decreases, showing that non-parallel effects become weaker for larger Reynolds numbers as expected. The missing terms in Eq. (19) considered by Bertolotti et al.20 compared to Eq. (18) are assumed to be of O(R−2), so the difference of the growth rate between Eqs. (19) and (18) is distinct when R = 998 but can be ignored when R is large.

TABLE I.

Comparison of eigenvalues obtained by different methods, where α = αr + i is the wavenumber obtained by solving the eigenvalue problem and α ̃ = α ̃ r + i α ̃ i is the corrected wavenumber defined in Eq. (23).

Method R ω αr αi × 102 α ̃ r α ̃ i × 1 0 2 α i / α ̃ i 1 ( % )
Jordinson24   998  0.112 200  0.308 600  −0.570 000       
LST  998  0.112 200  0.308 586  −0.572 708  0.308 586  −0.572 708   
EOS  998  0.112 200  0.307 408  −0.612 419  0.307 408  −0.612 419   
Bertolotti20   998  0.112 200  0.309 226  −1.151 846  0.307 297  −0.627 421  83.58 
EEV1  998  0.112 200  0.307 936  −1.246 716  0.306 814  −0.640 337  94.70 
EEV2  998  0.112 200  0.305 123  −1.493 983  0.307 232  −0.641 845  132.76 
LST  103  0.097 653  0.274 231  −0.740 209  0.274 231  −0.740 209 
Bertolotti  103  0.097 653  0.276 213  −0.964 522  0.272 754  −0.751 270  28.39 
EEV1  103  0.097 653  0.274 346  −1.008 233  0.272 691  −0.762 054  32.30 
EEV2  103  0.097 653  0.268 497  −1.197 989  0.272 739  −0.765 106  56.58 
LST  104  0.033 298  0.147 876  −1.220 273  0.147 876  −1.220 273   
Bertolotti  104  0.033 298  0.148 207  −1.366 281  0.147 686  −1.224 331  11.59 
EEV1  104  0.033 298  0.148 051  −1.366 283  0.147 682  −1.225 500  11.49 
EEV2  104  0.033 298  0.148 291  −1.445 845  0.147 685  −1.225 664  17.96 
LST  105  0.010 469  0.077 740  −0.812 477  0.077 740  −0.812 477   
Bertolotti  105  0.010 469  0.077 831  −0.851 569  0.077 720  −0.812 664  4.79 
EEV1  105  0.010 469  0.077 815  −0.851 535  0.077 720  −0.812 712  4.78 
EEV2  105  0.010 469  0.077 847  −0.878 676  0.077 722  −0.812 946  8.09 
Method R ω αr αi × 102 α ̃ r α ̃ i × 1 0 2 α i / α ̃ i 1 ( % )
Jordinson24   998  0.112 200  0.308 600  −0.570 000       
LST  998  0.112 200  0.308 586  −0.572 708  0.308 586  −0.572 708   
EOS  998  0.112 200  0.307 408  −0.612 419  0.307 408  −0.612 419   
Bertolotti20   998  0.112 200  0.309 226  −1.151 846  0.307 297  −0.627 421  83.58 
EEV1  998  0.112 200  0.307 936  −1.246 716  0.306 814  −0.640 337  94.70 
EEV2  998  0.112 200  0.305 123  −1.493 983  0.307 232  −0.641 845  132.76 
LST  103  0.097 653  0.274 231  −0.740 209  0.274 231  −0.740 209 
Bertolotti  103  0.097 653  0.276 213  −0.964 522  0.272 754  −0.751 270  28.39 
EEV1  103  0.097 653  0.274 346  −1.008 233  0.272 691  −0.762 054  32.30 
EEV2  103  0.097 653  0.268 497  −1.197 989  0.272 739  −0.765 106  56.58 
LST  104  0.033 298  0.147 876  −1.220 273  0.147 876  −1.220 273   
Bertolotti  104  0.033 298  0.148 207  −1.366 281  0.147 686  −1.224 331  11.59 
EEV1  104  0.033 298  0.148 051  −1.366 283  0.147 682  −1.225 500  11.49 
EEV2  104  0.033 298  0.148 291  −1.445 845  0.147 685  −1.225 664  17.96 
LST  105  0.010 469  0.077 740  −0.812 477  0.077 740  −0.812 477   
Bertolotti  105  0.010 469  0.077 831  −0.851 569  0.077 720  −0.812 664  4.79 
EEV1  105  0.010 469  0.077 815  −0.851 535  0.077 720  −0.812 712  4.78 
EEV2  105  0.010 469  0.077 847  −0.878 676  0.077 722  −0.812 946  8.09 

Figure 2 shows the distribution of the leading-order eigenvectors ( | u ˆ 0 | component) obtained by different methods for R = 998 and ω = 0.1122, and its distortion | u ˆ 1 | obtained by EEV1 is also given, whose amplitude is enlarged by a factor 102 for comparison purpose. As is illustrated, the leading-order eigenvectors obtained by LST, EOS, and EEV1 are almost identical to the graphical precision, and they reach the peak value at ym ≈ 0.44. The distortion | u ˆ 1 | in EOS and LST is zero, but is nonzero in EEV1 with its maximum value being about 0.6 × 10−2; the real part of u ˆ 1 , u ˆ 1 r , is about −0.6 × 10−2 at ym (while the imaginary part u ˆ 1 i is much smaller). However, such a small distortion of the eigenvector changes the growth rate substantially. The reason becomes clear if we note that the correction to the wavenumber due to distortion in Eq. (23) can be written as

α ˆ = α ˆ r + i α ˆ i = i 1 u ˆ u ˆ x | y m i u ˆ 1 u ˆ 0 | y m = u ˆ 1 i i u ˆ 1 r , where u ˆ 0 ( y m ) = 1 .
(26)

Since αi is rather small (of order 10−2), a small distortion u ˆ 1 r may change the growth rate by a amount comparable with the growth rate itself. That is why αi obtained by solving the eigenvalue problem is αi = − 0.012 467 16, but the corrected growth rate becomes α ̃ i = 0 . 006 403 37 . When R increases, the difference between αi and α i ̃ becomes smaller, as is displayed in the last column of Table I.

FIG. 2.

The eigenfunctions obtained by different methods for R = 998 and ω = 0.1122. (a) The leading-order eigenfunction | u ˆ 0 | and (b) the distortion of the shape function u ˆ 1 (which has been enlarged by a factor of 100) in EEV1 and the local solution proposed by Bertolotti.20 

FIG. 2.

The eigenfunctions obtained by different methods for R = 998 and ω = 0.1122. (a) The leading-order eigenfunction | u ˆ 0 | and (b) the distortion of the shape function u ˆ 1 (which has been enlarged by a factor of 100) in EEV1 and the local solution proposed by Bertolotti.20 

Close modal

Figure 3 shows the variation of the growth rate α ̃ i with the frequency F and the Reynolds number R. The predicted growth rates by EOS are somewhat higher than the LST result, and the unstable region is wider. Figure 3(a) indicates that at a given, R the growth rate obtained by EEV1 is smaller than that by LST for lower frequencies, but larger for higher frequencies, indicating that non-parallelism of the mean flow stabilizes T-S waves with lower frequencies while enhances T-S waves of higher frequencies.

FIG. 3.

Variation of the growth rate α ̃ i with (a) frequency F (R = 103) and (b) the Reynolds number R (F = 60). The differences of the growth rates predicted by EOS, Bertolotti,20 and EEV1 from that by LST: (c) variation with frequency F (R = 103) and (d) variation with the Reynolds number R (F = 60).

FIG. 3.

Variation of the growth rate α ̃ i with (a) frequency F (R = 103) and (b) the Reynolds number R (F = 60). The differences of the growth rates predicted by EOS, Bertolotti,20 and EEV1 from that by LST: (c) variation with frequency F (R = 103) and (d) variation with the Reynolds number R (F = 60).

Close modal

Similarly, Figure 3(b) shows that for a given frequency, the growth rate is reduced by non-parallel-flow effects before it reaches its peak value, and increases after that. The integral contribution of the growth rate to the amplitude would be changed slightly but the positions where the same amplitude is reached shifts downstream due to non-parallel-flow effects.

As was mentioned earlier, non-parallel-flow effects consist of two parts. In order to assess the relative role of each part, we display in Figures 3(c) and 3(d) the differences between the growth rates predicted by EOS and EEV1 and that by LST. Note that the label “EEV1(distorted)” refers to the prediction by EEV1 with the direct effect excluded. The curve labelled as “EOS-LST” represents the direct effect of non-parallelism in the mean flow. The difference is always positive, implying a stabilizing role. On the other hand, the curve labelled as “EEV1(distorted)-LST” measures the effect of shape distortion, which is found to stabilize/destabilize T-S waves with lower/higher frequencies or at lower/higher Reynolds numbers, respectively. Obviously, this effect is greater for most parameters. The curve labelled as “EEV1-LST” represents the sum of the two effects.

Figure 4 plots the neutral curves predicted by EOS and EEV1. They are compared with a set of earlier results that have been obtained by different methods, including the LST,24 the non-parallel theory,11 the residual minimization method,25 and DNS.15 As is illustrated, the result by EOS is similar to that by the residual minimization method;25 they both predict a larger unstable range than LST does. The result by EEV1 overlaps with those by DNS and the non-parallel theory.11 The frequencies for both the upper and lower branches of the neutral curve at a given Reynolds number are higher than the LST result. For the upper branch of the neutral cure, all of the methods give the same conclusion that non-parallelism of the mean flow enhances T-S waves. However, for the lower branch, different methods give different and even opposite conclusions. T-S waves are enhanced according to EOS and the residual minimization method but are inhibited according to EEV1 and the non-parallel theory.11 As was mentioned earlier, EOS is incomplete, and so its prediction should be discarded.

FIG. 4.

Comparison of the neutral curves obtained by EOS and EEV1 with those predicted by LST of Jordison (1970), and24 the non-parallel theory of Gaster (1974),11 the residual minimization method of Gaster (2000)25 and DNS of Fasel et al. (1990).15 

FIG. 4.

Comparison of the neutral curves obtained by EOS and EEV1 with those predicted by LST of Jordison (1970), and24 the non-parallel theory of Gaster (1974),11 the residual minimization method of Gaster (2000)25 and DNS of Fasel et al. (1990).15 

Close modal

The critical Reynolds numbers Rc predicted by LST and the non-parallel theory are 520 and 510, respectively. With the criterion of the inner maximum of |u′|, EEV1 gives Rc = 520 as LST does, but the corresponding neutral mode has a higher frequency. The critical Reynolds number is of little relevance to transition since T-S waves in this relatively high-frequency range undergo very limited accumulated amplification. Transition is caused instead by the growth of relatively low-frequency modes.

Characteristics of three-dimensional T-S waves have been studied in a series of experiments, which accumulated a comprehensive set of data.5,6,26,27 Guided by these measurements, we now perform LST, EEV1, and EEV2 calculations. Figure 5 displays the growth rates obtained by EEV1 and EEV2 for an oblique wave with F = 60 and β = 0.2 at different Reynolds numbers. For comparison, the corresponding results of LST and DNS are also presented. Considerable discrepancy exists between the LST and DNS results, indicating that non-parallelism has a rather strong effect on the three-dimensional wave. The first solution of EEV1 corresponding to the eigenvalue having the largest −αi is closer to that of DNS than LST does. The prediction of EEV2 coincides with that of DNS.

FIG. 5.

Growth rates α ̃ i obtained by EEV1, EEV2, LST, and DNS at different Reynolds numbers for F = 60 and β = 0.2 (normalized by the displacement thickness at R = 1000).

FIG. 5.

Growth rates α ̃ i obtained by EEV1, EEV2, LST, and DNS at different Reynolds numbers for F = 60 and β = 0.2 (normalized by the displacement thickness at R = 1000).

Close modal

In our DNS, we find that the inlet condition may have a significant effect on DNS results. Figure 6 shows the local growth rates computed by DNS for three different inlet disturbances, given by LST, EEV1, and EEV2, respectively. For a two-dimensional mode (Figure 6(a)), if LST is used as the inlet condition in the subcritical regime, a mild adjustment takes place. However, when EEV solutions are used, the adjustment disappears. For an oblique mode (Figure 6(b)), the growth rates collapse to a common asymptotic growth for R > 1150. However, when a LST normal mode is imposed (as is often done in DNS), considerable adjustment occurs in a region near the inlet, which is somewhat surprising despite that occurrence of adjustment is well known when the disturbance at the inlet is non-modal. Here, the adjustment occurs because without accounting for non-parallelism, the LST solution is not compatible with the non-parallel mean flow. With the EEV2 solution as the inlet condition, adjustment no long takes place, suggesting that EEV2 accurately approximates the asymptotic growth rate. The presence of adjustment phase sends a cautionary note: when DNS results with LST normal mode as inlet condition are employed to validate eigenvalue calculations, one should select the data sufficiently downstream the inlet, discarding those near the inlet. Such a care was exercised in the comparisons in the present paper.

FIG. 6.

Local growth rates obtained by DNS using different inlet conditions for (a) F = 60 and β = 0 and (b) F = 60 and β = 0.2, where β is the spanwise wavenumber normalized by the displacement thickness at R = 1000, the Reynolds number at the inlet.

FIG. 6.

Local growth rates obtained by DNS using different inlet conditions for (a) F = 60 and β = 0 and (b) F = 60 and β = 0.2, where β is the spanwise wavenumber normalized by the displacement thickness at R = 1000, the Reynolds number at the inlet.

Close modal

Characteristics of three-dimensional T-S waves are calculated using LST, EEV1, and EEV2 for the same parameters as in the experiment of Kachanov and Michalke.5 The results are displayed in Figure 7, where comparison is made with our DNS result as well as with the experimental data. Figures 7(a) and 7(b) show that for a wide range of wave propagation angles, all methods give virtually the same wavenumbers and phase velocities, which are in excellent agreement with the DNS and measurement, indicating that the wavenumber and phase velocity are insensitive to non-parallel-flow effects. However, non-parallelism strongly affects the growth rates of oblique waves, more so for larger wave angles, as is shown in Figure 7(c), in which the growth rates α ̃ i are calculated at the inner peak location of | u ˆ | . The predictions of EEV1 and EEV2 exhibit good agreement with the DNS. They are also closer to the experimental data than the LST prediction for modes with wave angle σ > 30°. However, the theoretical results are only in qualitative agreement with the experiment. In laboratory, it is more convenient to measure the disturbance amplitude and phase at a given distance from the wall. Figure 7(d) shows the growth rates α ̃ i ( y ) calculated at a given transverse location y = 1.273. A better agreement with experiments is observed for large wave angles.

FIG. 7.

Dependence of characteristics of oblique T-S waves on the wave angle σ (in degree) for R = 1129 and F = 94.1. (a) The (normalized) phase velocity cx(σ)/cx(0), (b) the (normalized) streamwise wavenumber α ̃ r ( σ ) / α ̃ r ( 0 ) , (c) the growth rate α ̃ i obtained at ym using | u ˆ 0 | max , and (d) the growth rate α ̃ i ( y ) obtained at a constant location y = 1.273. The experimental data are from Kachanov and Michalke.5 

FIG. 7.

Dependence of characteristics of oblique T-S waves on the wave angle σ (in degree) for R = 1129 and F = 94.1. (a) The (normalized) phase velocity cx(σ)/cx(0), (b) the (normalized) streamwise wavenumber α ̃ r ( σ ) / α ̃ r ( 0 ) , (c) the growth rate α ̃ i obtained at ym using | u ˆ 0 | max , and (d) the growth rate α ̃ i ( y ) obtained at a constant location y = 1.273. The experimental data are from Kachanov and Michalke.5 

Close modal

We have also calculated the local growth rates at different transverse locations for Reynolds number R = 1129, frequency F = 94.1, and different wavenumber β, corresponding to different wave angle σ. The variation with y is shown in Figure 8, where the shape function is also plotted in order to help interpret the results. In order to avoid overly crowded plots, we choose to present EEV2 results. While LST gives a growth rate independent of y, the growth rates obtained by EEV2 vary considerably with y in the region y < 4. In particular, in the vicinity of the local minimum of | u ˆ | , an extremely rapid variation occurs so that the local growth rate appears almost discontinuous. This near “discontinuity” arises possibly due to a nearly vanishing | u ˆ | , which “magnifies” 1 u ˆ u ˆ x . The location of the inner maximum of | u ˆ | is about y = 0.4 for planar waves, and shift upwards to y = 0.8 for fully oblique modes. For oblique modes, the growth rate α ̃ i ( y ) exhibits considerable variation near the inner maximum of | u ˆ | , and so the growth rate of oblique waves is more sensitive to the location y than that of planar waves, which makes accurate measurement of α ̃ i ( y m ) more difficult. The location of the outer maximum appears near y = 3, also most independent of σ. Interestingly, the growth rate calculated at this location using EEV2 (or EEV1) is the same as that given by LST for all wave angles. This feature is consistent with, and possibly explains, the experimental finding5 that the growth rates measured at the outer maximum agree better with the parallel-flow theory than those at the inner maximum do. Note that the predictions by EEV2 are in excellent agreement with our DNS results. The strong dependence of the local growth rate on the transverse position, especially the near “discontinuity,” was also observed in the experiment of Kachanov and Obolentseva.6 The experimental data are exacted from Figure 3 in their paper, and compared with our theoretical prediction in Figure 8. The error bars indicate that there is considerable uncertainty in the extracted data because Figure 3 contains many curves tangled together. Nevertheless, the overall agreement is rather good.

FIG. 8.

Variation of the local growth rate α ̃ i ( y ) with the transverse position y for R = 1129, F = 94.1 and different spanwise wavenumber β or wave angle σ: (a) σ = 0°, (b) σ = 21°, (c) σ = 38.9°, and (d) σ = 52.1°. The experimental data were taken from Kachanov and Obolentseva.6 

FIG. 8.

Variation of the local growth rate α ̃ i ( y ) with the transverse position y for R = 1129, F = 94.1 and different spanwise wavenumber β or wave angle σ: (a) σ = 0°, (b) σ = 21°, (c) σ = 38.9°, and (d) σ = 52.1°. The experimental data were taken from Kachanov and Obolentseva.6 

Close modal

In this paper, a systematic local and non-perturbative approach for studying non-parallel-flow effects on spatial instability of three-dimensional shear flows is presented using a three-dimensional boundary layer as paradigm. The streamwise variation of the mean flow affects the growth rate not only directly but also indirectly through the shape distortion of the perturbation. These effects are accounted for by expanding the mean-flow velocities as well as the perturbation shape as Taylor series. This leads to a sequence of extended eigenvalue problems, which are designated as EEVn with n referring to the number of terms retained in the Taylor series. These eigenvalue problems can be solved numerically, and the growth rate and wavelength of a disturbance with a given frequency can immediately be calculated. The approach is rather simple and efficient because although the extended eigenvalue problem EEVn involves, when discretized, a matrix whose dimension increases by n-fold compared with the usual eigenvalue problem for a parallel flow, the computational cost increases only moderately owing to the sparse structure of the matrix.

Unlike the perturbative WKB/multi-scale type of methods, which are valid only if non-parallelism appears as a small (higher-order) correction to the growth rate, the present approach is non-perturbative and the above restriction can be lifted. It is therefore applicable to flows where non-parallelism influences the leading-order growth rate of the perturbation. The approach is local and capable of predicting directly and accurately the asymptotic growth rate based solely on the local property of the mean flow. This is in contrast with the nonlocal PSE method, where the growth rate depends on the history of the mean flow and the disturbance, and undergoes adjustment unless the initial condition has already taken non-parallelism into account rather accurately. Furthermore, the present approach can, in principle, be extended to make the error arbitrarily small, while there exists a finite error in the PSE methodology because the step size for streamwise marching has to be greater than a cutoff value.

The extended eigenvalue problems, EEV1 and EEV2, were solved for the Blasius boundary layer. The growth rates obtained are found to be in excellent agreement with the prediction by the non-parallel theory of Gaster,11 with the DNS results of Fasel et al.,15 and our own. Comparison with experimental data was also made, and good qualitative agreement was observed. Though quantitative agreement is not perfect, the present approach (EEV1 and EEV2) performed appreciably better than other methods for oblique modes with propagation angles exceeding 30°.

For the Blasius boundary layer, our calculations reaffirm the fact that non-parallelism hardly affects the phase velocities and wavelengths of T-S modes, but significantly alters their growth rates, especially those of oblique modes. The direct effect of non-parallelism is found to be destabilizing, whereas the indirect effect through shape distortion turns out to be two-fold. At a fixed location (i.e., Reynolds number), it is stabilizing/destabilizing for low/high-frequency modes, and for a fixed frequency, it is stabilizing at relatively low Reynolds numbers but switches to a destabilizing role at high Reynolds numbers. Overall, the indirect effect through shape distortion is greater than the direct effect in the majority of the parameter space.

The growth rate depends strongly on the transverse location y for oblique T-S waves, including at the inner peak position of the streamwise velocity of the disturbance, suggesting that accurate positioning of the probe is needed in measurement if the data are to be compared with theoretical or numerical results. The growth rate calculated at the outer maximum is found, for all the parameter values considered, to coincide with the LST prediction (and hence could not reflect the effects of non-parallelism).

The approach has clear advantage over existing methods, and its further extension and application seem rather promising. As was mentioned earlier, it should be particularly suited for investigating cross-flow instability in three-dimensional boundary layers, where non-parallelism near the leading edge region affects the growth rate at leading order. It would also be interesting to extend the approach to supersonic boundary-layer instability, where lower-branch (two-dimensional) modes are unknown to be profoundly influenced by nonparallel-flow effects according to large-Reynolds-number asymptotic analysis (see, e.g., Ref. 28).

This work was supported by the National Natural Science Foundation (Nos. 11332007 and 11172204), the Tianjin Natural Science Foundation (No. 15JCYBJC19500), and the Open fund from State Key Laboratory of Aerodynamics (No. SKLA201401). The authors would like to thank Professor Jisheng Luo of Tianjin University for valuable discussions and Professor Yury Kachanov for providing the papers containing relevant experimental data.

When applied to a parallel mean flow, EEV1 (18) reduces to

M ϕ 0 ϕ 1 A ˆ p B ˆ p 0 A ˆ p ϕ 0 ϕ 1 = 0 0 .
(A1)

Let α be the eigenvalue of A ˆ p with the associated eigenvector ϕ, namely, A ˆ p ϕ = 0 , which is the O-S equation. Now, since | M ˆ | = | A ˆ p | 2 , M ˆ has the eigenvalue α with multiplicity two and a single associated eigenvector (ϕ, 0)T. For a non-parallel mean flow, the eigenvalue α of multiplicity two splits into two different eigenvalues, α1 and α2, of EEV1, similar to what was noted in Ref. 20. Obviously, for weakly non-parallel flows both α1 and α2 are close to α.

In general, when EEVn formulation is applied to a parallel mean flow, the operator will have the same eigenvalue α but with multiplicity (n + 1) and a single associated eigenfunction (ϕ, 0, …, 0)T. For a non-parallel flow, the eigenvalue α will split into (n + 1) different perturbed eigenvalues α0, …, αn in the equations EEVn with different eigenvectors.20 However, the physical solution of the eigenvalue problem ought to be unique, namely, the corrected wavenumber α ̃ defined in Eq. (23) should be the same. Next, we shall provide computational evidence to confirm that this is indeed the case.

Figure 9 shows the growth rates obtained by EEV1 and EEV2 for planar and oblique T-S waves. For comparison purpose, the results of LST and DNS are also presented. As is illustrated, EEV1 and EEV2 have two and three eigenvalues, respectively, which are appreciably different. In the case of planar modes, the corrected (i.e., physical) growth rates α ̃ i obtained by EEV1 are indistinguishable (Figure 9(a)), and EEV2 also yields the same physical growth rate from three distinct eigenvalues (Figure 9(b)). For planar T-S waves, EEV1 gives sufficiently accurate prediction in excellent agreement with the DNS prediction, and EEV2 is not needed. In the case of oblique modes, the two corrected growth rates α ̃ i obtained by EEV1 are much closer to DNS as well as to each other. However, there remain appreciable differences among them (see Figure 9(c)), implying that the EEV1 prediction is not yet accurate enough for oblique modes. Proceeding to EEV2, we find that all three eigenvalues lead to practically identical predictions in good agreement with the DNS result, as is shown in Figure 9(d). In practice, it suffices to choose the eigenvalue α with the largest (−αi). Our experience with the calculations suggests that this leads to a somewhat more accurate result compared with using other eigenvalues.

FIG. 9.

The corrected growth rates α ̃ i obtained from different eigenvalues −αi using EEV1 and EEV2 for planar and oblique modes with F = 60. (a) β = 0.0, EEV1, (b) β = 0.0, EEV2, (c) β = 0.2, EEV1, (d) β = 0.2, EEV2. The spanwise wavenumber β is normalized by the displacement thickness at R = 1000.

FIG. 9.

The corrected growth rates α ̃ i obtained from different eigenvalues −αi using EEV1 and EEV2 for planar and oblique modes with F = 60. (a) β = 0.0, EEV1, (b) β = 0.0, EEV2, (c) β = 0.2, EEV1, (d) β = 0.2, EEV2. The spanwise wavenumber β is normalized by the displacement thickness at R = 1000.

Close modal

We also calculated damped discrete modes as well as continuous spectra, and the results are displayed in Figure 10. For each of three discrete modes (marked by red symbols), EEV1 and EEV2 yield two and three eigenvalues, respectively, which are clustered around that predicted by LST. For the unstable mode, the two eigenvalues of EEV1 converge to a common value after correction (23). For each of two damped modes, the two eigenvalues of EEV1 become much closer. Extending the calculation to EEV2, all three eigenvalues collapse to a single value for the first damped mode. Unfortunately, EEV2 was unable to provide reliable eigenvalues for the second damped mode possibly because the latter was too close to the continuous spectra, which reside along the almost vertical lines. Similar to its discrete counterparts, each continuous mode splits into two in EEV1 formulation, but they overlap after correction (23). These continuous modes are noted to have phase speeds close to unity and were referred to as vortical modes. They have traditionally been regarded as representing vortical disturbances in the free stream. This long-held view was found to be a misconception by the recent work of Wu and Dong,29 who showed that these modes are non-physical.

FIG. 10.

Spectra calculated by LST, EEV1 and EEV2 for R = 1000. (a) Eigenvalues of LST, EEV1 and EEV2 without correction. (b) Eigenvalues of EEV1 and EEV2 with the correction given by Eq. (23).

FIG. 10.

Spectra calculated by LST, EEV1 and EEV2 for R = 1000. (a) Eigenvalues of LST, EEV1 and EEV2 without correction. (b) Eigenvalues of EEV1 and EEV2 with the correction given by Eq. (23).

Close modal
1.
H. L.
Reed
,
W. S.
Saric
, and
D.
Arnal
, “
Linear stability theory applied to boundary layers
,”
Annu. Rev. Fluid Mech.
28
(
1
),
389
428
(
1996
).
2.
E.
Reshotko
, “
Boundary-layer stability and transition
,”
Annu. Rev. Fluid Mech.
8
(
1
),
311
349
(
1976
).
3.
G.
Schubauer
and
H.
Skramstad
, “
Laminar boundary-layer oscillations and transition on a flat plate
,”
J. Res. Natl. Bur. Stand.
38
,
251
(
1947
).
4.
J. A.
Ross
,
F. H.
Barnes
,
J. G.
Burns
, and
M. A. S.
Ross
, “
The flat plate boundary layer. Part 3. Comparison of theory with experiment
,”
J. Fluid Mech.
43
(
04
),
819
832
(
1970
).
5.
Y. S.
Kachanov
and
A.
Michalke
, “
Three-dimensinoal instability of flat-plate bounday layers: Theory and experiment
,”
Eur. J. Mech. B-Fluid.
13
(
4
),
401
422
(
1994
).
6.
Y. S.
Kachanov
and
T. G.
Obolentseva
, “
Development of three-dimensional disturbances in the Blasius boundary layer 3. Nonparallism effects
,”
Thermophys. Aeromech.
5
(
3
),
331
338
(
1998
).
7.
T.
Herbert
, “
Parabolized stability equations
,”
Annu. Rev. Fluid Mech.
29
(
1
),
245
283
(
1997
).
8.
M. D. J.
Barry
and
M. A. S.
Ross
, “
The flat plate boundary layer. Part 2. The effect of increasing thickness on stability
,”
J. Fluid Mech.
43
(
04
),
813
818
(
1970
).
9.
M.
Bouthier
, “
Stabilite lineaire des ecoulements presque paralleles
,”
J. Mec.
11
,
599
621
(
1972
).
10.
C. H.
Ling
and
W. C.
Reynolds
, “
Non-parallel flow corrections for the stability of shear flows
,”
J. Fluid Mech.
59
,
571
591
(
1973
).
11.
M.
Gaster
, “
On the effects of boundary-layer growth on flow stability
,”
J. Fluid Mech.
66
(
03
),
465
480
(
1974
).
12.
W. S.
Saric
and
A. H.
Nayfeh
, “
Nonparallel stability of boundary layer flows
,”
Phys. Fluids
18
,
945
950
(
1975
).
13.
Th. L.
Stijn
and
A. I.
Vooren
, “
On the stability of almost parallel boundary layer flows
,”
Comput. Fluids
10
(
4
),
223
241
(
1982
).
14.
J. B.
Thomas
and
J. M.
Philip
, “
Boundary layer stability calculations
,”
Phys. Fluids
30
(
11
),
3351
3358
(
1987
).
15.
H.
Fasel
and
U.
Konzelmann
, “
Non-parallel stability of a flat-plate boundary layer using the complete Navier–Stokes equations
,”
J. Fluid Mech.
221
,
311
347
(
1990
).
16.
D. G.
Crighton
and
M.
Gaster
, “
Stability of slowly diverging jet flow
,”
J. Fluid Mech.
77
,
397
413
(
1976
).
17.
C. C.
Lin
, “
On the stability of two-dimensional parallel flows
,”
Q. J. Mech. Appl. Maths
3
,
117
142
(
1946
).
18.
F. T.
Smith
, “
On the non-parallel flow stability of the Blasius boundary layer
,”
Proc. R. Soc. London, Ser. A
366
(
1724
),
91
109
(
1979
).
19.
R. J.
Bodonyi
and
F. T.
Smith
, “
The upper branch stability of the Blasius boundary layer, including non-parallel flow effects
,”
Proc. R. Soc. London, Ser. A
375
(
1760
),
65
92
(
1981
).
20.
F. P.
Bertolotti
,
Th.
Herbert
, and
P. R.
Spalart
, “
Linear and nonlinear stability of the Blasius boundary layer
,”
J. Fluid Mech.
242
,
441
474
(
1992
).
21.
M. R.
Turner
and
P. W.
Hammerton
, “
Asymptotic receptivity analysis and the parabolized stability equation: A combined approach to boundary layer transition
,”
J. Fluid Mech.
562
,
355
381
(
2006
).
22.
Z. F.
Huang
,
H.
Zhou
, and
J. S.
Luo
, “
Direct numerical simulation of a supersonic turbulent boundary layer on a flat plate and its analysis
,”
Sci. China, Ser. G
48
,
626
640
(
2005
).
23.
Z. F.
Huang
,
W.
Cao
, and
H.
Zhou
, “
The mechanism of breakdown in laminar-turbulent transition of a supersonic boundary layer on a flat plate-temporal mode
,”
Sci. China, Ser. G
48
,
614
625
(
2005
).
24.
R.
Jordinson
, “
The flat plate boundary layer. Part 1. Numerical integration of the OrrC-Sommerfeld equation
,”
J. Fluid Mech.
43
(
04
),
801
811
(
1970
).
25.
M.
Gaster
, “
On the growth of waves in boundary layers: A non-parallel correction
,”
J. Fluid Mech.
424
,
367
377
(
2000
).
26.
Y. S.
Kachanov
and
T. G.
Obolentseva
, “
Development of three-dimensinoal disturbances in the Blasius bounday layers 1. Wave-trains
,”
Thermophys. Aeromech.
3
(
3
),
225
243
(
1996
).
27.
Y. S.
Kachanov
and
T. G.
Obolentseva
, “
Development of three-dimensinoal disturbances in the Blasius bounday layers 2. Stability characteristics
,”
Thermophys. Aeromech.
4
(
4
),
373
384
(
1997
).
28.
F. T.
Smith
, “
On the first-mode instability in subsonic, supersonic or hypersonic boundary layers
,”
J. Fluid Mech.
198
,
127
153
(
1989
).
29.
X.
Wu
and
M.
Dong
, “
On continuous spectra of the Orr-Sommerfeld/Squire equations and entrainment of free-stream vortical disturbances
,”
J. Fluid Mech.
732
,
616
659
(
2013
).
30.
A.
Butler
and
X.
Wu
, “
Non-parallel-flow effects on stationary crossflow vortices at their genesis
,”
in Proc. IUTAM Symposium on Laminar-Turbulent Transition, Rio de Jeneiro, Brazil,
edited by
M.
Medeiros
(
Elsevier
,
2015
, in press).