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.

## I. INTRODUCTION

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/2}*x*^{∗}, much smaller than *x*^{∗} when the Reynolds number $Re= U \u221e \u2217 x \u2217 / \nu \u221e \u2217 $ is large, where $ U \u221e \u2217 $ is the free-stream velocity and $ \nu \u221e \u2217 $ 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 \u0304 $ 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 Re^{1/2} > 600,^{3,4} but appreciable discrepancy exists for Re^{1/2} < 600 due to non-parallel-flow effects.^{1} It is recognized that non-parallelism is more significant for oblique T-S waves^{5,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, $\u2202 u \u0304 /\u2202x$ and $ v \u0304 $ and (b) the influence of the streamwise distortion of the local eigenfunction. The contribution (a) can easily be accounted for by including $\u2202 u \u0304 /\u2202x$ and $ v \u0304 $ 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. Gaster^{11} 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 Gaster^{11} was remarkably accurate (see below). Saric and Nayfeh^{12} and Crighton and Gaster^{16} 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 Re^{1/8} times the boundary-layer thickness, and the well-known triple-deck structure emerges.^{17,18} Using this formalism, Smith^{18} 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 Bouthier^{9} and Gaster^{11} 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 Smith^{19} 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 *x* − *x*_{0}, for example, the disturbance shape function *ϕ*(*x*, *y*) = *ϕ*_{0}(*y*) + (*x* − *x*_{0}) *ϕ*_{1}(*y*) + *O*((*x* − *x*_{0})^{2}), as was done in Ling and Reynolds,^{10} but the main differences are that *x*_{0} 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*[(*x* − *x*_{0})], 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 Hammerton^{21} 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 Konzelmann^{15} 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 Gaster^{11} 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 *x _{a}* with

*x*is then approximated by a Taylor expansion,

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 Konzelmann^{15} 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.

## II. FORMULATION

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 \u221e \u2217 , \rho \u221e \u2217 $, and $ \mu \u221e \u2217 $ 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 $ \delta 0 \u2217 $ and $ \delta 0 \u2217 / U \u221e \u2217 $, respectively, where $ \delta 0 \u2217 $ is the displacement thickness of the boundary layer at a typical location $ x 0 \u2217 $; the resulting dimensionless coordinates and time will be denoted as (*x*, *y*, *z*) and *t*, respectively. The velocity and pressure are normalized by $ U \u221e \u2217 $ and $ \rho \u221e \u2217 U \u221e \u2217 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

The dimensionless three-dimensional incompressible N-S equations read

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.

### A. Perturbations and linear instability

Let the three-dimensional mean flow field be denoted by $ Q \u0304 ( x , y ) = ( u \u0304 , v \u0304 , w \u0304 , p \u0304 ) $. 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

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

here we have put

### B. Eigenvalue problem

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 \u0304 $ is also present). The perturbation **Q**′(*x*, *y*, *z*, *t*) can be written as

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

where

here we have put

The matrices $ A \u02c6 p $, $ B \u02c6 p $, and $ C \u02c6 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 \u02c6 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.

### C. The O-S equation

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

which holds at an arbitrary location *x _{a}* 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

*α*(

*x*) and the corresponding eigenvector

_{a}*ϕ*(

*x*,

_{a}*y*) at a location

*x*for a given

_{a}*β*and

*ω*. The eigenvector

*ϕ*(

*x*,

_{a}*y*) is usually normalized by $ u \u02c6 ( x a , y m ) $, where

*y*is the location at which $| u \u02c6 |$ 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.

_{m}### D. Extended eigenvalue problems

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 *x _{a}* is approximated by a Taylor expansion, namely,

in which the streamwise variation of *α* has been absorbed into the *O*((*x* − *x _{a}*)

^{2}) term in the expansion of the shape function. It follows that ∂

*α*/∂

*x*in Eq. (11) can be set to zero since

*α*=

*α*(

*x*) is a constant in the vicinity of

_{a}*x*. This term is also of smaller size than other terms representing non-parallelism when the Reynolds number is large.

_{a}Similarly, the variation of the mean flow is accounted for by expanding its velocities, $ u \u0304 $, $ v \u0304 $, and $ w \u0304 $, as Taylor series. Correspondingly, the operator expands as

and similar expansions are performed for $ B \u02c6 p $ and $ A \u02c6 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 Reynolds^{10} 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

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,

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}*ϕ*/∂*x*^{2} 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* = *x _{a}*.

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

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

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}*ϕ*/∂*x*^{2}. Similar to LST and EEV1, this eigenvalue problem can be solved directly with the local information of the mean flow at *x* = *x _{a}*.

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

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 \u02c6 , v \u02c6 , w \u02c6 , p \u02c6 ) \u2192 ( 0 , 0 , 0 , 0 ) $ as *y* → −∞.

### E. Corrected wavenumber and growth rate

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 *x _{a}*, 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,

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 \u2032 = q \u02c6 exp { \u222b x \alpha d x} $, once

*α*and

*q*′ are found by solving the extended eigenvalue problem, the corrected wavenumber is given by

where $ \alpha \u02c6 $ stands for the correction to the wavenumber due to the distortion of the shape function $ q \u02c6 ( y ) $. Usually, the streamwise velocity $ u \u02c6 $ of the perturbation is chosen and *y* is taken to be the peak position of $| u \u02c6 |$ since $ u \u02c6 $ is the quantity that can be measured rather accurately. Unless stated otherwise, this is how the corrected wavenumber $ \alpha \u0303 $ 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

In addition to its growth rate $\u2212 \alpha \u0303 i $, other characteristics of interest are its propagation angle *σ* and the streamwise phase velocity *c _{x}*, which are defined, respectively, as

## III. NUMERICAL METHODS

### A. The eigenvalue problems

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 \u0304 ( 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* = *x _{a}* 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 \u02c6 $ at the wall is, as usual, remedied by discretizing the normal pressure gradient, $\u2202 p \u02c6 /\u2202y$, 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

where *k*_{1} and *k*_{2} are parameters. The default grid parameters in Eq. (25) are selected as *k*_{1} = 100, *k*_{2} = 60, the position of the outer boundary *y*_{max} = 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 $ \alpha \u0303 $, 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 −

*α*. A further discussion of multiplicity and uniqueness is given in the appendix.

_{i}### B. Direct numerical simulations

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 \u221e \u2217 =34m/s$, $ \rho \u221e \u2217 =1.225kg/ m 3 $, $ p \u221e \u2217 =101\u2009325Pa$, $ \mu \u221e \u2217 =1.79\xd71 0 \u2212 5 Pa\u2009s$, and $ c \u221e \u2217 =340m/s$, which give the unit Reynolds number *Re*^{∗} = 1.1634 × 10^{6}/m and a Mach number *M* = 0.1. For a given Reynolds number $R= Re \u2217 \delta 0 \u2217 $, there is a corresponding displacement thickness $ \delta 0 \u2217 $ and the streamwise $ x 0 \u2217 $, 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 *A*_{0} = 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 *x _{i}* =

*x*

_{0}+

*i*Δ

*x*(0 ≤

*i*≤

*I*), where Δ

*x*= (

*x*−

_{e}*x*

_{0})/

*I*with

*x*

_{0}and

*x*denoting the positions of the inlet and outlet of the computation domain, respectively, and

_{e}*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

*C*= 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.

_{fl}## IV. RESULTS

### A. Verification

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 Konzelmann^{15} and the theoretical result of Gaster;^{11} the latter itself has been confirmed by experiments,^{3,4} DNS^{15} 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.

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.

### B. Eigenvalues and eigenvectors

Table I displays the results obtained by different methods, in which the result of LST by Jordinson^{24} and the local solution by Bertolotti^{20} is listed as a reference. When *R* = 998, the growth rate $\u2212 \alpha \u0303 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.

Method . | R
. | ω
. | α
. _{r} | α × 10_{i}^{2}
. | $ \alpha \u0303 r $ . | $ \alpha \u0303 i \xd71 0 2 $ . | $ \alpha i / \alpha \u0303 i \u22121 ( % ) $ . |
---|---|---|---|---|---|---|---|

Jordinson^{24} | 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 | |

Bertolotti^{20} | 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 | 10^{3} | 0.097 653 | 0.274 231 | −0.740 209 | 0.274 231 | −0.740 209 | |

Bertolotti | 10^{3} | 0.097 653 | 0.276 213 | −0.964 522 | 0.272 754 | −0.751 270 | 28.39 |

EEV1 | 10^{3} | 0.097 653 | 0.274 346 | −1.008 233 | 0.272 691 | −0.762 054 | 32.30 |

EEV2 | 10^{3} | 0.097 653 | 0.268 497 | −1.197 989 | 0.272 739 | −0.765 106 | 56.58 |

LST | 10^{4} | 0.033 298 | 0.147 876 | −1.220 273 | 0.147 876 | −1.220 273 | |

Bertolotti | 10^{4} | 0.033 298 | 0.148 207 | −1.366 281 | 0.147 686 | −1.224 331 | 11.59 |

EEV1 | 10^{4} | 0.033 298 | 0.148 051 | −1.366 283 | 0.147 682 | −1.225 500 | 11.49 |

EEV2 | 10^{4} | 0.033 298 | 0.148 291 | −1.445 845 | 0.147 685 | −1.225 664 | 17.96 |

LST | 10^{5} | 0.010 469 | 0.077 740 | −0.812 477 | 0.077 740 | −0.812 477 | |

Bertolotti | 10^{5} | 0.010 469 | 0.077 831 | −0.851 569 | 0.077 720 | −0.812 664 | 4.79 |

EEV1 | 10^{5} | 0.010 469 | 0.077 815 | −0.851 535 | 0.077 720 | −0.812 712 | 4.78 |

EEV2 | 10^{5} | 0.010 469 | 0.077 847 | −0.878 676 | 0.077 722 | −0.812 946 | 8.09 |

Method . | R
. | ω
. | α
. _{r} | α × 10_{i}^{2}
. | $ \alpha \u0303 r $ . | $ \alpha \u0303 i \xd71 0 2 $ . | $ \alpha i / \alpha \u0303 i \u22121 ( % ) $ . |
---|---|---|---|---|---|---|---|

Jordinson^{24} | 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 | |

Bertolotti^{20} | 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 | 10^{3} | 0.097 653 | 0.274 231 | −0.740 209 | 0.274 231 | −0.740 209 | |

Bertolotti | 10^{3} | 0.097 653 | 0.276 213 | −0.964 522 | 0.272 754 | −0.751 270 | 28.39 |

EEV1 | 10^{3} | 0.097 653 | 0.274 346 | −1.008 233 | 0.272 691 | −0.762 054 | 32.30 |

EEV2 | 10^{3} | 0.097 653 | 0.268 497 | −1.197 989 | 0.272 739 | −0.765 106 | 56.58 |

LST | 10^{4} | 0.033 298 | 0.147 876 | −1.220 273 | 0.147 876 | −1.220 273 | |

Bertolotti | 10^{4} | 0.033 298 | 0.148 207 | −1.366 281 | 0.147 686 | −1.224 331 | 11.59 |

EEV1 | 10^{4} | 0.033 298 | 0.148 051 | −1.366 283 | 0.147 682 | −1.225 500 | 11.49 |

EEV2 | 10^{4} | 0.033 298 | 0.148 291 | −1.445 845 | 0.147 685 | −1.225 664 | 17.96 |

LST | 10^{5} | 0.010 469 | 0.077 740 | −0.812 477 | 0.077 740 | −0.812 477 | |

Bertolotti | 10^{5} | 0.010 469 | 0.077 831 | −0.851 569 | 0.077 720 | −0.812 664 | 4.79 |

EEV1 | 10^{5} | 0.010 469 | 0.077 815 | −0.851 535 | 0.077 720 | −0.812 712 | 4.78 |

EEV2 | 10^{5} | 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 \u02c6 0 |$ component) obtained by different methods for *R* = 998 and *ω* = 0.1122, and its distortion $| u \u02c6 1 |$ obtained by EEV1 is also given, whose amplitude is enlarged by a factor 10^{2} 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 *y _{m}* ≈ 0.44. The distortion $| u \u02c6 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 \u02c6 1 $, $ u \u02c6 1 r $, is about −0.6 × 10

^{−2}at

*y*(while the imaginary part $ u \u02c6 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

_{m} Since *α _{i}* is rather small (of order 10

^{−2}), a small distortion $\u2212 u \u02c6 1 r $ may change the growth rate by a amount comparable with the growth rate itself. That is why

*α*obtained by solving the eigenvalue problem is

_{i}*α*= − 0.012 467 16, but the corrected growth rate becomes $ \alpha \u0303 i =\u22120.006\u2009403\u200937$. When

_{i}*R*increases, the difference between

*α*and $ \alpha i \u0303 $ becomes smaller, as is displayed in the last column of Table I.

_{i}### C. Growth rates

Figure 3 shows the variation of the growth rate $\u2212 \alpha \u0303 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.

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.

### D. Neutral curve

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.

The critical Reynolds numbers *R _{c}* predicted by LST and the non-parallel theory are 520 and 510, respectively. With the criterion of the inner maximum of |

*u*′|, EEV1 gives

*R*= 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.

_{c}### E. Three-dimensional T-S waves

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.

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.

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 $\u2212 \alpha \u0303 i $ are calculated at the inner peak location of $| u \u02c6 |$. 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 $\u2212 \alpha \u0303 i ( y ) $ calculated at a given transverse location *y* = 1.273. A better agreement with experiments is observed for large wave angles.

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 \u02c6 |$, 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 \u02c6 |$, which “magnifies” $ 1 u \u02c6 \u2202 u \u02c6 \u2202 x $. The location of the inner maximum of $| u \u02c6 |$ 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 $\u2212 \alpha \u0303 i ( y ) $ exhibits considerable variation near the inner maximum of $| u \u02c6 |$, 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 $\u2212 \alpha \u0303 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 finding^{5} 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.

## V. SUMMARY AND CONCLUSIONS

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

## Acknowledgments

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.

### APPENDIX: MULTIPLICITY OF EIGENVALUES

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

Let *α* be the eigenvalue of $ A \u02c6 p $ with the associated eigenvector *ϕ*, namely, $ A \u02c6 p \varphi =0$, which is the O-S equation. Now, since $| M \u02c6 |=| A \u02c6 p | 2 $, $ M \u02c6 $ 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 EEV

*n*with different eigenvectors.

^{20}However, the physical solution of the eigenvalue problem ought to be unique, namely, the corrected wavenumber $ \alpha \u0303 $ 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 $\u2212 \alpha \u0303 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 $\u2212 \alpha \u0303 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.

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.