Spectral analysis of parabolized stability equations

Computers & Fluids

Vol. 26. No. 3. pp. 219-291. 1997 0 1997 Ekvier Science Ltd Printed in Great Britain. All rights reserved 004s7930/97 $17.00 + 0.00

Pergamon

PII: so045-793q%)ooo44-8

SPE’CTRAL

ANALYSIS

OF PARABOLIZED EQUATIONS

STABILITY

FE1 LI and MUJEEB R. MALIK* High Technology Corporation,

Hampton, VA 23666, U.S.A.

(Received 20 May 1996; in revised form

9 September 1996)

Abstract-In recent years, the parabolized stability equations (PSE) approach has gained popularity for transition studies in boundary-layer flows. In this paper, the physical origin of ill-posedness of PSE in the primitive-variable form is studied by spectral analysis. In subsonic flows, a continuous spectrum representing the upstream-propagating acoustic waves is shown to be responsible for the ill-posedness of PSE. In supersonic flows, acoustic waves no longer contribute to ill-posedness despite the existence of a subsonic region within the boundary layer. In this case, ill-posedness is caused by some discrete modes with upstream influence. Semi-discretized models of PSE are similarly studied. It is shown that successful implementalion of PSE to physical problems can only be achieved through manipulation of the spectra of the PSE operator. A generalized condition for the stable marching solution is devised. 0 1997 Elsevier Science Ltd. All rights reserved.

1. INTRODUCTION

In boundary-layer type of flows, various instability mechanisms may cause transition to turbulence. These instabilities often take the form of waves with distinct wavelengths and frequencies. Furthermore, if small perturbations to the laminar flow field are triggered by a source (oscillating or otherwise) at some upstream location, the development of these perturbations downstream seems to be determined by what happens at the source. This suggests that an initial-value problem in space may be a good model of the instability mechanism, i.e. initial perturbations are specified at some streamwise location in the boundary layer and their evolution is traced by marching downstream. One such model is the parabolized stability equations (PSE) for slowly-varying, nonparallel boundary layers. PSE was first used by Herbert [l] in streamfunction form (see also Bertolotti et al. [I!]).The primitive-variable form of PSE was used by Bertolotti and Herbert [3] and Chang et al. [4] for compressible flows. In practice, the PSE method proves to be very efficient and successful and is used to solve many difficult problems. However, it was noticed by Chang et al. [4] that the solution procedure fails for the primitive-variable form of PSE if the discretization step-size in the streamwise direction becomes sufficiently small. This type of solution failure is the consequence of ill-posedness of the problem, which was also observed in solving parabolized Navier-Stokes equations (PNS) [5]. The physics of this phenomenon will be analyzed in this paper. For simplicity, we will first consider a two-dimensional boundary layer. The extension of the analysis to three-dimensional flows is straightforward and will be considered later. Let x and y be the streamwise and wall-normal coordinates, respectively, and U and V are the associated mean velocity components. We impose a small disturbance of frequency w on the mean flow and denote it by cp = (u,v,p)’ which represent the velocities and pressure in the usual notation. We can write cp in the following way: 1

cp(x,y,t) = $(x,y)exp

(S i

a,(x’)dx’ - iot vll

+ C.C.

(1)

where $ = (ir,Q)‘. Most of the streamwise oscillations and growth of cp are absorbed by the complex wavemmber a, (x), leaving +(x,y) a slowly-varying function of x. Therefore, terms involving 8*I(//dx* can be neglected. The two-dimensional primitive-variable form of PSE

*Corresponding

author. 279

F. Li and M. R. Malik

280

equations [6] can be written as

aii

ai,

where C = irx,U - io + c$lRe - idu,/dx and K = U - 2icr,lRe. The Reynolds number, Re, is based on some length scale proportional to the boundary-layer thickness at some fixed streamwise station. The condition that determines the choice of c(, can take several forms, for example

s

x.. u*-aii

0

ax

dy = 0

(5)

where ic* is the complex conjugate of i2. Initial conditions are specified at some upstream location. Equation (5) is discussed in more detail by Herbert [7] and Malik et al. [8]. The well-posedness of a problem requires that (i) the solution exists, (ii) is unique, and (iii) it must continuously depend on the initial and boundary conditions, i.e. the change of the solution must be small for small perturbations to the initial or boundary conditions. The usual way of classifying partial differential equations with real coefficients is the characteristic analysis. An equation is said to be elliptic if some characteristics are complex. However, the characteristic method does not, in general, address the issue of continuous dependence on initial conditions. Analyzing the continuous dependence of solutions on the initial conditions is nontrivial for a system of equations with variable coefficients and, therefore, such a task is often carried out by freezing the coefficients of the equations. However, as Kreiss and Lorenz [9] pointed out, the well-posedness properties of the frozen equations may not necessarily represent those of the ‘unfrozen’ ones. Furthermore, for equations with complex coefficients, some complex characteristics may not represent ellipticity and care must be taken in interpreting these characteristics. The terms ‘ellipticity’ and ‘elliptic effects’ are repeatedly used in this paper. We will, therefore, provide a working definition of these terms in the context of the present analysis. Whenever an admissible solution of the governing equations is associated with upstream propagation of information, we will say that this solution represents ellipticity. Haj-Hariri [lo] and Herbert [7] applied the characteristic method to PSE by freezing the coefficients of the equations. They found three complex characteristics, two of which are associated with acoustic wave propagation and the other is believed to be viscosity related. These were all interpreted as representing ellipticity. Their methods of analysis essentially consist of replacing X/ax, &jay and a’iijay’ in the equations by @, Ai2and A2ir,respectively, and finding p as a function of I, i.e. p =f(n). If this function is complex, the equation system is said to be elliptic. This works well for equations with real coefficients, but may run into trouble when the equations have complex coefficients. For example, consider the following equation with complex constant coefficients:

d@Y)

=fti>

(7)

wheref(Y) is a well-behaved function whose Fourier transform exists. The unique solution to this problem is given by $(%Y) =

(8)

For x < 0, the real part of the exponent in the above integral becomes positive and a small, but

Parabolized stability equations

281

general, perturbation to f(y) will cause the integral to ‘blow up’. However, in the region x > 0, a small perturbation to_&) will change the solution by a small amount. Therefore, the solution indicates that the information specified at x = 0 propagates towards the region where x > 0. It can also be shown that the above solution is unique. Thus, all three requirements of well-posedness are satisfied. This system is in fact parabolic according to the definition of Kreiss and Lorenz [9]. However, by applying the characteristic method just described, one would find a complex characteristic p = (1 + i) A2and conclude, erroneously, that the problem modeled by equations (6) and (7) is elliptic and cannot be solved as an initial-value problem in x. Li and Malik [1 1, 121carried out well-posedness analysis of PSE by using the method of Kreiss and Lorenz [9]. The coefficients of the equations are frozen. This method essentially consists of replacing Xi/ax, &?/ay and a’b/ay* in the equations by pit, Ui’i and - ,12ir,respectively, and finding p as a function of A, i.e. p =f(n). If RealCf(A)) < 0 as ];l]+co, the problem is well-posed. Otherwise, unbounded solutions exist and a small perturbation to the initial conditions may result in very large changes of the solution at any finite x. Application of this method to equation (6) yields p = - (1 + i)A2 and, hence, the problem is well-posed. Li and Malik [1 1, 121 found three unbounded solutions for PSE, which are related to the three characteristics of Haj-Hariri [lo] and Herbert [7]. Two of the solutions were acoustic waves and the other was found to represent ellipticity only in a small neighborhood of the wall. We will show later that the third solution does not cause ellipticity if the analysis is carried out without the frozen-coefficient assumption. The frozen-coelhcient methods are used in the analyses discussed. However, the rapid variation of the mean flow in the wall-normal direction casts some doubt on the validity of the frozen-coefficient assumption. More sophisticated methods are necessary for unveiling the physical mechanisms that cause ellipticity. In this paper, we will study the PSE without freezing the coefficients in the ,wall-normal direction. The physical origins of ellipticity are individually identified and ways of managing or suppressing ellipticity are analyzed. Two other conclusions of the previous analyses are that the upstream propagation of information is not related to the convective terms in the governing equations [lo, 71 and that the sign of the streamwise velocity in the boundary layer does not affect the mathematical nature of PSE [7]. These claims will also b’e examined in this paper. Finally, we explain what is meant by ‘parabolization’. In their book, Kreiss and Lorenz [9] provide the following definition for a parabolic system: consider an equation of the form

where x and y ar’e real; 4 is a complex vector of dimension n; A, B and C are complex matrices of dimension n x n. This system is defined as parabolic if all n complex eigenvalues of matrix A have positive real parts. In general, PSE cannot be written in the form of equation (9). Analysis in this paper shows that the PSE system is still elliptic. Hence, the term ‘parabolization’ in the present context only refers to the process of reducing the elliptic effect to a minimum by modifying the linear Navier-Stokes operator, while still retaining the dominant physical character of the desired disturbance mode. The plan of th.e paper is as follows: in Section 2, we deal with the incompressible two and three-dimensional boundary layers. In Section 3, the theory is extended to compressible flows. Section 4 contains an example used to verify the theory for compressible flows. The paper is concluded in Section 5.

2. INCOMPRESSIBLE BOUNDARY LAYERS We will carry out spectral analysis of the two-dimensional linearized Navier-Stokes equations (LNS) with three goals in mind: (i) introduce the methodology of spectral analysis and hence, pave the way for detailed analysis of PSE; (ii) discuss the discrete and continuous spectra for the LNS operator under thie parallel flow assumption; and (iii) identify the subset of the spectra which gives rise to elliptic effect, and hence, facilitate comparisons with the spectra of PSE later in this section and clarify what exactly is meant by ‘reducing’ or ‘eliminating’ ellipticity. Many of the results derived for incompressible equations can be easily extended to compressible equations.

F. Li and M. R. Malik

282

2.1. Linearized incompressible Navier-Stokes

equations

Here, we consider linearized Navier-Stokes equations for parallel flows. The analysis of nonparallel flows is beyond the scope of the present work. Bertolotti [13] analyzed wave reflections due to streamwise nonparallelism using model equations. We start by transforming the time derivatives in the governing equations. In a systematic analysis of an initial-value problem in time, the generalized Laplace-Fourier transform must be carried out. However, in the spatial problem modeled by PSE, no such formalism is explicitly used. Hence, in our linearized Navier-Stokes equations, we shall simply ‘factor out’ the time derivatives by using the function e-‘“’ with w denoting the fixed frequency. The following equations are thus obtained for parallel flows:

-

icou+Ug+d(iv= -iov+

ap

-~+$v’u

dy

(10)

U+jf = - -ap -I- $V2u

(11)

ay

2+&=O ay

(12)

Suppose we desire to trace the development of the least stable mode in the boundary layer starting at some location x = x0. Initial conditions need to be specified at such a location. We assume here that the initial conditions consist of the least stable mode and some noise (perturbations). The boundary conditions are as follows: u(x,O) = v(x,O) = 0

(13)

lim u(x,y) = lim v(x,y) = 0 1-a J-JO

(14)

and

It is well-known that the above initial-value problem cannot be solved because of ‘elliptic effects’. In fact, Salwen and Grosch [14] attempted to solve this problem and encountered difficulties associated with upstream-propagating disturbances. Here, in order to illustrate the difficulties, we will consider this initial-value formulation. For an initial-value problem, it is natural to use the one-sided Laplace transform in x

s cc

(&$,j)’ =

(u,v,p)Te

ia-dx

(15)

0

Thus, equations (lo)--(12) become the following set of nonhomogeneous equations

ordinary

differential

(17) d5 . ,. KXu+ dy = U” where C = i(crU - w) + a’/Re and uo, vo, po, &,/ax, &,/ax are the imposed initial conditions at x = x0. The boundary conditions are that ir and B vanish at the wall and that they are bounded as y-+00. Without going into details of how the above ODES can be solved, we merely assume that the solution, (ii$,j)‘, is at our disposal. Our remaining task is to obtain the solution (~,v,p)~ by inverse Laplace transform. Hence, we need to find the inversion contour which, in the complex a-plane should run below all the poles and branch cuts of the solution. Therefore, our first step must be to locate these poles and branch cuts in the complex a-plane. These poles and branch cuts are given

Parabolized stability equations

283

by the spatial eigenspectra of equations (16~( 18) with zero right-hand side. Henceforth, we shall use the term ‘poles’ interchangeably with ‘discrete eigenvalues’, and the term ‘branch cuts’ interchangeably with ‘continuous spectra’. For incompressible boundary-layer flow, there exist a finite number of discrete eigenvalues and a finite number of continuous spectra in the complex a-plane. In a properly posed analysis (e.g. [ 15, 161)using generalized Laplace-Fourier transform for an initial-value problem in time, the eigenspectra can be correctly interpreted. When a point disturbance located at some streamwise location, say x,, creates a perturbation field in the flow, the poles and branch cuts sufficiently above the real a-axis contribute to the disturbance field downstream of x, and those sufficiently below the real a-axis contribute to that upstream. In the neighborhood of the real a-axis, the upstream or downstream association of the poles and branch cuts can also be found, but the method is complicated and since it is described in Refs [15] and [16], we will not repeat it here. A clear illustration of the upstream or downstream association of the poles can be found in Fig. 5 of Ashpis and Reshotko [l:!]. We now denote the particular downstream-associated pole with the most negative imaginary part in the complex a-plane by D,, and call it the least stable or the most unstable mode. Our task is therefore to trace the development of D, downstream. However, we must be aware that, in our problem, the specification of the initial conditions at x = x0 implicitly assumes that a small disturbance at so:me location, x,, only influences the flow field downstream of x,, but not upstream, hence, all the poles and branch cuts in the complex a-plane are interpreted by the differential operator in equations (lOk(12) as downstream modes. Consequently, should there exist upstream-associated poles or branch cuts below D, in the complex a-plane, the operator would wrongly interpret them as downstream growing modes with growth rates greater than that of the D, mode and any perturbation to the initial condition at x,, having components in these modes would contaminate the correct solution. We will call these modes the ‘parasite’ modes. As we will see very soon, these parasite modes do exist in the present problem. Suppose we ha.ve identified all the finite poles (including D, in particular) and we need now to find all the branch cuts or the continuous spectra. Gustavasson [17] was the first to show that the branch cut in the Laplace inversion plane introduces the continuous spectrum for the temporal case and subsequent ,works have determined the continuous spectra in the spatial case in the same manner [1416]. 1t will become apparent that the continuous spectra play an important role in our subsequent analysis, therefore, we will describe, in some detail, how these continuous spectra are obtained. For a parallel boundary layer, the velocity becomes constant (U, = 1) at large distances from the wall. Therefore, the homogeneous solutions in the freestream can be obtained analytically. The general solution for y+ cc takes the form

(19)

where A,,*= + ,,& L and A,, = f iRe(a - co) + a’. The solution for y+cc must not be infinite and hence 2, = A, = 0. In order to avoid double-valued functions we must have Real(s) > 0 and Real( iRe(a - co) + a’) > 0. This leads to the requirement that a2 and iRe(a - co) + a2 must not be any negative real number. Hence, the equations governing the location of the branch cuts in the complex a-plane are a= a=

-

iRe + iRed

*iv + 4(q2 - ioRe)/Re2 2

(20)

(21)

where 0 < q < a). These represent four branch cuts in the complex a-plane, which are plotted schematically in IFig. 1, together with the discrete eigenvalue D,. The first and second branch cuts (20) have branch points at the origin (0,O) and extend along the imaginary axis to positive and negative infinities,, respectively. The third branch cut has its branch point at approximately (w,O + )

F. Li and M. R. Malik

284

and extends to + ice and the fourth branch cut has its branch point at approximately ( - w, - Re) and extends to - ioo. Both these latter branch cuts asymptotically tend to the imaginary a-axis at large distances from the real axis. We denote these branch cuts by C,, Cz, C, and C,, respectively. The correct physical interpretations of these four branch cuts are as follows: C, and C, both have zero wave number, i.e. Real(a) = 0, and hence have infinite phase speed (c = o/Real(a)), therefore, they represent acoustic wave propagation in the incompressible flow limit. C, is associated with downstream propagation of small disturbances with phase speed equal to 1 (the freestream speed) and hence represents slow decaying convective disturbances. C, is associated with rapidly decaying, upstream viscous diffusion. However, as explained before, C, and C, are wrongly interpreted as downstream growing solutions in the above initial-value problem and, therefore, a Laplace inversion contour cannot be found and hence our problem is ill-posed. It might be tempting to argue that, if the initial conditions do not have any components in C, and C, we may still trace the correct solution. However, this is virtually impossible because any random perturbation to the initial conditions having components in these two branch cuts will cause the solution procedure to fail. Furthermore, since both C, and C,, extend to - ice in the a-plane, there should be no bounds on the growth rate of the parasite solutions associated with these branch cuts. Hence, our problem is ill-posed because the solution does not continuously depend on the initial conditions. A simple example of ill-posedness of this kind can be found in the Hadamard problem when the Laplace equation is solved as an initial-value problem ([ 181,p. 90). Since C, and C, are the main culprits for ellipticity, it may be possible that ellipticity can be eliminated or, at least, reduced by altering the differential operator so that C, or C, no longer appear in the spectra. However, we must not change the physical character of the mode, D,, which we wish to trace. Parabolized stability equations partially achieve this goal. We should note that the solutions given by equation (19) do not necessarily vanish as y+ 00. In physical space, the solutions are expressed in terms of integrals of equation (19) along the branch cuts. The vanishing boundary conditions as y -+ co given by (14) are satisfied through cancellation. This can be illustrated by considering the inverse Fourier transform of some function, f(y), that vanishes at ]y]-+co

None of the base functions, eiLy,vanishes as Jy]+co, but the integral does. 2.2. Spectral analysis of incompressible PSE We have already introduced the parabolized stability equations in Section 1. Now we will carry out the spectra1 analysis of PSE. In order to simplify the problem, we again analyze a parallel boundary layer and, in addition, we assume that u,--the complex wavenumber of mode D,-is fixed. We should keep in mind that PSE solves for the slowly-varying shape function of the instability mode that is traced. In a parallel flow, such a shape function does not vary at all. We shall see what happens if the most unstable mode, D,, is specified as initial condition at x = x0 and traced downstream. The PSE is obtained through not only the substitution given by equation (l), but also an approximation, i.e. neglecting the second derivatives with respect to x. The differential operator has been altered. The new equations for a parallel flow are obtained by setting V(x,y) = au/ax = 0 in equations (2)-(4) with initial conditions consisting of mode D, and some noise, and the boundary conditions are as before. Once again, we will apply the Laplace transform to these equations and investigate whether an inversion contour can be found. Here we use 5 as the Laplace variable to avoid confusion with a used in the previous case. The resulting equations are (C + iSK)ic + $C (C + i(K

+ i(a, + 03 + $

-

&

j& 3

$

= Ku, + po

= KU0

(22) (23)

Parabolized stability equations

285

,. i(a, + 5); + ’ = 240 dy

where C = i(a,U .- o) + af/Re and K = U - 2ia,/Re. As before, we must find all the poles and branch cuts of the solution to equations (22X24), which are given by the eigenspectra of these equations with vanishing right-hand side. In particular, if we set 5 = 0, (22)-(24) become identical to (16x18) for a = a,. Since the pole D, is represented by a, in the complex a-plane, it follows that l = 0 is a pole in the complex r-plane. Furthermore, the eigenfunctions in the two cases are exactly the same. Therefore, the substitution given by equation (1) has shifted the pole D, in the a-plane to the origin of the l-plane, and the dropping of the second derivatives with respect to x did not have any effect on this mode in the parallel limit, We will continue to denote this mode by D,. The wavenumber of mode D,, Real([), is zero, implying that the shape function does not oscillate downstream. Furthermore, the growth rate of the shape function, - Imag(& is also zero, i.e. the amplitude of the shape function stays constant. The oscillation and growth of the physical wave are absorbed in the complex wavenumber, a,. The other poles in the complex t-plane can be computed by solving equations (22k(24) with zero right-hand side as an eigenvalue problem. The poles near the origin of the c-plane ([ N O(1) or smaller) can actually be obtained approximately by simply shifting the corresponding poles in the a-plane by a,. This can be explained as follows. If the neglected second derivatives with respect to x were added ‘back to the equations, then (22) and (23) would each contain an extra term of order, t2/Re. The eigenvalues would be the same as those for the linear Navier-Stokes equations shifted by a,. Hence, deleting the second derivatives with respect to x creates a perturbation of order <‘/Re to the diIferentia1 operator. While small but random perturbations may produce large changes in the eigenvalues of the hydrodynamic stability operators [ 191, systematically deleting terms like [‘/Re will lead only to negligible changes for eigenvalues of order 1 or smaller [20]. Thus, these poles in the a-plane near the origin are approximately shifted by a, into the t-plane. The poles or branch cuts far away from the origin are modified, some even disappear. Let us not worry about the poles for the moment since we know, from analyzing the linear Navier-Stokes equations, that the continuous spectra give rise to the most dangerous parasite solutions. The continuous spectra of PSE can be found in a similar manner as before. The branch cuts C, and C, in the complex a-plane are shifted by an amount of exactly a, to arrive at their locations in the complex <-plane and they still extend to ice and - ice, respectively. The branch cut C, is shifted approximately by a, and extends to ioo and the branch cut C, has disappeared. We shall still denote these branch cuts by C,, C, and C,. The exact equations for these remaining branch cuts are 5=

-a,+&

Fig. 1. Schematic plot of the four branch cuts in the incompressible linear Navier-Stokes spectra.

(25)

Fig. 2. Schematic plot of the three branch cuts in the incompressible PSE spectra.

F. Li and M. R. Malik

286

(= t=

-a,-iiv]

i(a: + q’) + (0 - a,)Re Re(U, - 2iaJRe)

(26) (27)

where 0 < ye< co. We are using U, in the last expression instead of its value, 1, to facilitate our discussion. These branch cuts are schematically shown in Fig. 2. Since C, runs into the lower half of the complex a-plane and extends to - ice, we immediately realize that it will prevent the inversion of the Laplace transform, causing ellipticity. The branch cuts C, and C, extend into the upper half of the complex g-plane and, therefore, do not cause problems for the inversion of the Laplace transform. The disappearance of branch cut C, is due to the fact that the differential operator has been modified by neglecting the second derivatives with respect to x and, as a result, the streamwise diffusion of the parasite solution no longer occurs. Thus, we may say that the ellipticity has been ‘reduced’. However, the parasite solution associated with branch cut C, has unbounded growth and will cause our solution procedure to fail. Hence, the remaining ellipticity is by no means weak. From equation (26) for branch cut C,, it may appear that the upstream propagation of information is not related to the convective terms in PSE because it can be verified that M, and iv in the expression do not come from streamwise convective terms in equations (22)-(24), as stated by Haj-Hariri [lo] and Herbert 171. However, this is only true for acoustic waves in the incompressible limit. For compressible boundary layers, the convective velocity can be comparable to the speed of sound and, hence, convective terms in PSE will affect upstream-propagating acoustic waves. It can be shown that the branch cuts are closely related to the characteristics of PSE. As q + co, equations (25)-(27) for branch cuts C,, C, and C, take on similar forms to the three complex characteristics found by Haj-Hariri [lo] and Herbert [7]. In particular, the characteristic corresponding to C, was thought to be an upstream-associated viscous mode because the Reynolds number appears in the equation of the complex characteristic. In the well-posedness analysis of Li and Malik [l 1, 121 this mode was found to represent ellipticity only in the immediate neighborhood of the wall. According to the present analysis, C, does not represent ellipticity at all. It turns out that this mode is actually a downstream-associated convective mode. Analogous to equation (6) discussed in the Section 1, a complex characteristic of PSE, whose coefficients are complex, may not necessarily give rise to upstream influence. The characteristic corresponding to C, is complex, but it does not represent ellipticity. The analysis of Li and Malik [l 1, 121 was not based on characteristics, but suffers from the limitations of the frozen-coefficient assumption. Because of such an assumption, the freestream velocity U, in equation (27) was wrongly interpreted as the local velocity U(V) anywhere in the boundary layer, and hence, if Imag(a,) < 0, i.e. for a growing wave, and U(y) < - Imag(a,)/Re, e.g. very close to the wall, the sign of the imaginary part in equation (27) was reversed and C, was thought to cause ellipticity. It is also interesting to note that C, would extend into the lower half of the Laplace inversion plane and cause ellipticity if the sign of U, were reversed. Hence, we can only solve the problem in the direction of freestream convection, but not against it. The mathematical nature of branch cut C, is best illustrated by the following parabolic equation:

a4 -=

ax

1 a: + (W - a,)Re 824 Re(U, - 2iaJRe) zdy ‘jRe(U,2iaJRe) 4

where 4 is some model function. The only branch cut of equation (28) in the Laplace inversion plane is the same as that given by equation (27). It is clear that C, presents parabolicity and is the branch cut to which PSE owes its name. In practice, one effective way to ‘reduce’ the remaining ellipticity associated with C2 is to drop the pressure gradient term @/ax from equation (2) (see Refs [4] and [lo]). We multiply the term @/ax in equation (2) by a constant 0. If D = 1 the PSE is not modified. If cr = 0, we have two branch cuts remaining. We can also carefully let rr approach zero from above in a continuous fashion and examine the migration of the branch cuts. In this process, the branch cut C, is moved up in the a-plane and disappears at ( - co, + 00) branch cut C, is unchanged (27) and branch cut C, approaches a limit and its location in the t-plane is given by

287

Parabolized stability equations

t=

-~(l +j$)-k,(l - &)

where 0 < rl < co. The branch point of Cz is unaffected by the limiting process and remains at ( - %r, - IX,,).If the mode D, is unstable, i.e. cl,, < 0, equation (29) implies that branch cut C, extends into the lower half of the complex g-plane diagonally. The growth rate of the parasite solution associated with C, is still unbounded. Therefore, by setting @/ax = 0, we have eliminated a harmless branch cut C,, while the targeted branch cut, Cz, has been modified but not removed. However, we can see from (29) that the location of C, is moved farther away from D, and later analysis will show that the effect of this parasite solution is reduced in a semi-discretized system. 2.3. Semi-discretrzed

incompressible

PSE

In practice, the discretized PSE system is quite well-behaved if the marching step-size in x is sufficiently large, but still small enough to resolve the slowly-varying shape function of D,. However, we should not become complacent just because solutions can be obtained despite the ill-posedness of the original equations. We need to know how much physics is retained or lost by using a discretized model with large step-sizes and what happens to the parasite solutions. To analyze the phenomenon of stabilization by large step-sizes, we only need to discretize the PSE in the x-direction. The semi-discretized system is obtained by one-step backward difference with step-size 6x and the resulting equations are given below (30)

Gin,,

1 d*fi,+,

din+,

+ -

dy

K,

--7=Xvn Re dy

(32) where G = C + k:/hx, K = (U - 2iaJRe) and C = i(a,U - w) + afIRe. Boundary conditions are such that the velocities vanish at the wall and infinity. For such systems, the Laplace transform analysis used in the preceding cases is no longer applicable. Hence, we will use the Z-transform [21]. For the present purpose, we need only to introduce the definition of the transform and its inversion formula. For n = 0,1,2 ,..., etc., let {cp,} be a complex sequence. Assuming that there exists a positive real number A such that Iv,,+ ,/(pJ < A for n -*co, then, for some complex number z with ]z] > A, the following sum exists:

W) =

2 (Pzn

(33)

This is the Z-transform which maps the complex sequence {qn> to a function of the complex z-plane. The inversion formula is given below 1 (Pn= G

z”- ‘$(z)dz

(34)

The function I(i(z) is analytic for ]z] > A and its analytic continuation exists everywhere inside the circle ]zJ = A, except at the poles or on the branch cuts. Therefore, the inverse Z-transform can be expressed in terms of these poles and branch cuts. Furthermore, if a pole or part of a branch cut is outside the unit circle ]z] = 1, the associated solutions are growing ones. A useful property of the Z-transform is that, if I(/(z) is the Z-transform of the sequence {cp.}, then, z$(z) - zcp, is the Z-transform of the sequence {cp,+ ,}.

F. Li and M. R. Malik

288

We now apply the Z-transform

to equations (30)-(32) and the resulting equations are

(C + SKyi +

$6 +i(ct, +

e)j -

t

$ =

r,

(35)

1 d2i? di (C + iOK)u + j-j - jjg dy2 = r2 1

. i(cq + tl)ic + * = r, dy

(37)

where C = i(u,U - w) + uf/Re and K = U - 2iu,/Re. The right-hand side of equations (35)-(37), r,, r2, and r3 are functions of the initial conditions at x = x0, independent of z. The variable 8 is given by

If we set z = 1 in the complex z-plane, i.e. 0 = 0, we again recover equations (16x18) for a = a,. Hence, z = (1 ,O) is a pole in the complex z-plane, corresponding to mode D,. Since lz] = 1 for mode D,, this mode neither grows nor decays. Furthermore, this mode is independent of the discretization step-size. These properties are results of using a parallel flow model in which the shape function does not change with x. For weakly nonparallel flows, the change of the shape function is minimized (5) and the effect of large step-size on D, is small. Hence, the physics of D, is mostly preserved. The branch cuts C,, C, and C, in the z-plane are given, respectively, as 1 z = (1 - f&$x + r/6x) + &,6x

(39) (40)

and (Re - 2ia,) ’ = (Re - 2icq) + (v + af/Re + icq - io)bx

(41)

where 0 < ye< co. It can be easily verified that C, and C, fall inside the circle of unit radius, centered at the origin of the complex z-plane (for a,, < 0), i.e. they represent exponentially decaying solutions. Branch cut C, is bounded by a circle of radius l/la,dx] for any finite step-size 6x. If we choose the inversion contour to be a circle of sufficiently large radius so that all the poles and branch cuts are enclosed, the sum given by equation (33) is analytic outside the circle, i.e. the Z-transform is invertible, our problem becomes well-posed in the conventional sense. However, the solution on part of the branch cut grows exponentially (the part outside the unit circle in Fig. 3) for sufficiently small step-size 6x and will eventually contaminate the desired solution if the projection of the initial conditions to C2 is nonzero (however small). Therefore, for PSE to be a useful tool, it must satisfy stricter conditions than the conventional definition of well-posedness. Hence, an additional condition must be that C,, as well as the poles, reside on or within the unit circle. If la,,Sxl > 1 (Fig. 3), then the entire branch cut C, retreats completely into the unit circle while some poles may still be outside. A necessary condition for the minimum step-size is, therefore (42)

This is the same condition obtained by Li and Malik [ll, 121. We will discuss shortly the circumstances under which equation (42) is also a sufficient condition. This condition is derived by using one-step backward difference. If multi-step backward differences were used, the minimum step-size increases. However, one-step difference is very special in that the error is proportional to the product of the step-size and the second derivative, i.e. (r?‘/ax’)sx. If one-step difference were to cause large error for reasonable step-sizes, it must mean that the second derivative is large-a violation of the PSE assumption. One must then carefully examine the validity of solutions

Parabolized

stability

equations

289

Fig. 3. Schematic plot of the z-plane of the semi-discretized PSE. Solid lines with symbols situations where 6x < AX,,,,,. Dashed lines with symbols represent 6x > 6x,,..

represent

obtained in such cases. On the other hand, if one-step difference causes negligible changes in the solution for a range of sufficiently large step-sizes, it must imply that the PSE assumption holds. We note that C, and C, form a closed circle which prevents the integration contour from getting into the interior of the region bounded by this circle (Fig. 3). This difficulty can be overcome by the method of Ashpis and Reshotko [15] in which the branch points of C, and C, are separated by a small distance 6 so that the integrand is analytic in-between. After the integration contour has gone through, 6 is allowed to approach zero. Once the poles, and branch cuts are manipulated into the unit circle (with the exception of D, which is always Ion the circle), the desired solution can be traced. It is important to notice that all the unwanted solutions decay exponentially. The semi-discrete system (30t(32) does not admit algebraic solutions except for constant solutions (e.g. 0,). We further notice that (35x37) are identical to (22k(24), with 5 replaced by 13.Therefore, the entire complex r-plane is mapped onto the complex z-plane by the transformation

or, equivalently

We have already shown that mode D, is mapped from 5 = (0,O) to z = (1,0) and have found the locations of the c80rresponding continuous spectra in the c-plane and z-plane. Any poles or branch cuts that are outside the unit circle in the complex z-plane will represent parasite solutions. From (44), it is clear that this unit circle corresponds to a circle of radius 1/6x and centered at

Fig. 4. Correspondence

between

the two complex

planes.

290

F. Li and M. R. Malik A

c-plane

Fig. 5. Upstream-associated

modes must be outside of this circle for well-posedness.

5 = (0, - 1/6x) in the complex {-plane (see Fig. 4). We denote these circles by Cr and Cc, respectively. A point outside C, is mapped to a point inside C,, More importantly, a point inside Cs is mapped to a point outside C_. Therefore, parasite solutions will contaminate mode D,. Hence, the necessary and sufficient condition for marching stability is that C, must be sufficiently small (i.e. 6x must be sufficiently large) so that it contains no poles or branch cuts. In the case of the Blasius boundary layer, no poles exist inside C, once the acoustic branch cut C, is excluded by setting 6x > l/lol,J. The condition given by equation (42) is, therefore, also sufficient for marching stability. A more general result which takes into consideration both poles and branch cuts can be obtained as follows. Let B be a pole other than D, or a point on a branch cut whose real and imaginary parts in the c-plane are given by tUrand ~,i, respectively. The minimum step-size is reached when B rests on the circle C, in the <-plane. Figure 5 depicts this situation for rUr< 0 and tti < 0. From elementary geometry, we obtain 53, +

(

l2 5, + z

>

1 = 6x2

Hence (46) Since B can be any pole or any point on a branch cut, the minimum step-size is given by 6X,i,

=

-

max (

ttr

25ui +

5ti

>

(47)

where the maximum is taken over all poles and all points on the branch cuts. The condition given by equation (42) derives from this result if (26) is substituted for (c,,, t,i). When poles and branch cuts do not exist in the lower half of the complex c-plane i.e. 5.i > 0 for all poles and branch cuts, there is no step-size restriction. Let B be a pole (other than 0,) or a point on a branch cut with coordinates 5, = (cti, 5.i). We will define a nonlinear distance of the point B to D, (the origin in t-plane) by (48) Thus, it is in this sense that we measure the closeness of a pole to D,. In a semi-discretized system, the poles or branch cuts closest to D, are the most dangerous ones. If an upstream-associated pole

Parabolized stability equations

291

is very close to II;),,the step-size required will be very large. Furthermore, the collision of B and D, gives rise to absolute instability. Under these circumstances, the PSE method would fail. We have shown earlier that deleting +/ax in the PSE moves the branch cut C, away from D,. The above results explain why ellipticity due to C, is reduced in that case. If two-step backward difference is used in the marching direction, the corresponding transformation between the t-plane and z-plane is given by < = i 3z2 - 42 + 1 2z26x This transformation is not one-to-one. A point in the t-plane is mapped to two points in the z-plane. In particular, the origin of the l-plane (mode 0,) is mapped to z = (1,O) and z = (l/3,0); both are poles in the z-plane. Therefore, the initial conditions must be specified on the former and any noise on the latter will decay as (l/3)“. We will not go into any more details here other than stating that, for a. stable marching solution, it is sufficient to exclude upstream poles and branch cuts from a circle centered at 5 = (0, - 2/6x) with radius 2.3/6x. This is a slightly conservative condition. The step-size required for well-posedness for two-step difference is, therefore, larger than that for one-step difference. The analysis we have carried out is not limited to incompressible flows. In fact, for any linear partial differential equation system which is first order in x and whose coefficients are constant with respect to x, the poles and branch cuts in the two complex planes are related by equation (44) when one-step difference is used. Therefore, the minimum step-size required to eliminate parasite solutions (if they exist) is given by (47). This would spare us of deriving everything again for compressible flows. 2.4. Three-dimensional incompressible J¶OWS The extension of the theory to three-dimensional flows is straightforward. In the analysis we will not only allow the waves to be oblique, but also introduce an additional mean flow component, W = W(y), the spanwise velocity. This change in the mean flow has no effect on the acoustic branch cuts in the incompressible limit. The branch cut C2 in the r-plane is given by [=

-ci-iJm

(49)

where 0 < q < cc and /? is the spanwise real wavenumber (assumed positive without loss of generality). The branch point is at a finite distance below the real axis. The mode D, is at the origin. In a semi-discretized system, it is easy to verify that, if I/I + C(iilI Ic(J, the step-size restriction is still given by (42). Otherwise

Y/J+ 4il dxmm = (p + f&J2+ cx:,< &J In many practical situations where IpI >>la,,], this can be approximated

3. COMPRESSIBLE

BOUNDARY

by

LAYERS

The method of analysis for compressible boundary layers is the same as for incompressible ones and makes use of the spectra of the operators of equations for compressible flows. 3.1. Linearized compressible Navier-Stokes

equations

In the incompressible flow limit, the branch cut associated with upstream propagation of acoustic waves is mainly responsible for ill-posedness of the PSE formulation. For compressible flows with subsonic freestream Mach number most of the branch cut C, shown in Fig. 2 is continuously pushed to the left of the complex a-plane as the Mach number increases. In this process, the distance between C2 and D, becomes increasingly large and the numerical problems caused by C2 diminish. For M > 1, C, appears in the right half of the a-plane above the real axis. Thus, C, becomes a downstream-associated branch cut and no longer contributes to ill-posedness. Any

F. Li and M. R. Malik

292

remaining ellipticity must, therefore, come from poles or branch cuts that are left in the lower half of the cl-plane. Consequently, ill-posedness of an initial-value problem in a supersonic boundary layer takes on a totally different physical character from that in an incompressible flow. The eigenspectra of a parallel supersonic boundary layer were discussed in Balakumar and Malik [16], and the detailed derivation will not be repeated here. After applying the Laplace transform in x to the compressible linear Navier-Stokes equations, we are again faced with the task of finding an inversion contour in the complex a-plane. In addition to poles, there exist three viscosity-associated branch cuts in the lower half of the a-plane which extend to - ice. The imaginary parts of the branch points have very large negative values (O(Re)). These branch cuts prevent us from finding an inversion contour and hence represent ill-posedness. The locations of these branch cuts in the cl-plane are independent of the form of the boundary-layer profile near the wall, therefore, the initial-value problem would still be ill-posed even if the entire flow field were supersonic. The elliptic effects associated with these three branch cuts are not due to the existence of a subsonic layer near the wall, they are there because of viscosity. In a semi-discretized system, if the step-size is small enough (6x N O(l/Re)) to resolve these continuous modes, the initial-value problem will fail. In practice, these poles do not cause much trouble since they are far away from D,. There are also poles in the lower half of the cc-plane, which will be discussed later in the context of PSE. In the upper half a-plane, there are four branch cuts, two of which extend to + ice almost vertically. They do not interfere with the inverse Laplace transform. We denote these two branch cuts by C, and C5, respectively (schematically shown in Fig. 6). The former is the same as C, in the incompressible case discussed before, and C, is a new branch cut associated with entropy waves. The other two branch cuts are associated with acoustic wave propagations and reduce to their respective counterparts in the incompressible flow as the freestream Mach number approaches zero. We still denote them by C, and C,, respectively (schematically shown in Fig. 6). Both are now downstream modes. Therefore, acoustic waves associated with these branch cuts will not propagate upstream when the freestream is supersonic despite the existence of a subsonic region within the boundary layer. At this point, we can see the limitations of the frozen-coefficient method in analyzing well-posedness discussed in Section 1. When the coefficients of the equations are frozen in y, the analysis is done locally. For local Mach numbers less than unity, one would conclude that acoustic waves travel upstream and give rise to elliptic effects. It is important to remember, however, that the familiar wave equation governing the properties of acoustic waves is derived by neglecting

a-plane

Fig. 6. Schematic plot of the four downstream-associated branch cuts in the linear Navier-Stokes spectra for supersonic boundary layers.

Parabolized

stability

equations

293

viscosity and mean velocity gradient. Hence, we have

(

AP - w’p

- 2io 58P

>

=(I

-M3g

azp azp + ay’

where u is the frequency. This equation shows that sound waves can travel upstream when M c 1. In a boundary layer, neither viscosity nor mean velocity gradient is negligible. Therefore, equation (52) no longer represents the underlying physics of wave propagation inside the boundary layer. Results from analyses based on the frozen-coefficient assumption may indicate that some disturbances have upstream influence, but will not, in general, reveal the physical nature of such disturbances. The upstream influence inside the boundary layer is mainly due to upstream-associated discrete modes (poles). Whether some of these discrete modes can be called acoustic modes may be a matter of definition, but they certainly do not satisfy the above acoustic equation (52), not even approximately. We denote the mode we desire to trace by D, and the upstream decaying mode closest to D, by B,. We will show that B, is responsible for ill-posedness of PSE. In what follows, we will discuss the parasite solutions which represent ill-posedness in compressible flows. 3.2. Linearized compressible PSE The primitive-variable formulation of PSE for compressible flows is given by Chang et al. [4]. We will not present the details of the equations. We will analyze a quasi-three-dimensional boundary layer with edge Mach numbers M in the marching direction (x-direction) and M, in the spanwise direction, and apply the Laplace transform in x. The spectra of the compressible PSE operator can be obtained in the manner as described before. The Navier-Stokes spectra in the neighborhood of the mode D, in the a-plane are approximately shifted by a, into the PSE spectra in the t-plane (ie. for 4 N O(1) or smaller), where a, is the complex wave number of mode D, in the Navier-Stokes spectra. Of particular interest is how the acoustic branch cuts change as A4 is increased from 0. The determination of acoustic branch cuts in the compressible case is nontrivial. However, since these branch cuts are associated with the mean flow outside the boundary-layer edge, the inviscid theory results in a good approximation (except for A4 z 1). Thus, an approximate equation for the acoustic branch cuts in the t-plane is given by

r= -a,+

-COM2+/?MM,f& 1 - M2

(53)

where 0 < q < a), Q = (oM - j?Mr)’ - (1 - M2)(q2 + 8’) and fi is the spanwise wavenumber of the oblique wave. If M, M, and /l are set to 0, equations (25) and (26) for the incompressible two-dimensional wave are recovered. When the parameter q is increased from zero, the two branch cuts represented ‘by equation (53) initially run parallel to the real axis in the negative &-direction. When q2 > - /I’ + (wM - /3M,)2/(1 - M’), the expression under the square root, Q, becomes negative and the branch cuts turn sharply through 90”. One runs vertically into the upper half of the complex plane and the other into the lower half, which causes ellipticity. If viscous effects were included, the turning of the branch cuts from the real axis would happen smoothly. Therefore, according to equation (47)

hri,” =

1

(54)

WM2 - /IMM, % + 1 - M2 for semi-discretizled PSE until A4 is such that an upstream discrete mode becomes closer to D, in the sense of equation (48). For sufficiently large jI, the minimum step-size is actually smaller than that given by (54), similar to the incompressible case expressed by (50). However, it will not be presented here. For supersonic freestreams (M > I), the quantity Q is always positive, and the acoustic branch cuts run along the real axis of the complex plane. If viscous effects were included, the two branch cuts would gradually move into the upper half of the complex plane as shown

294

F. Li and M. R. Malik

c-plane

B, 0

Fig. 7. Schematic plot of the four downstream-associated branch cuts in the linear PSE spectra for supersonic boundary layers.

schematically by curves C, and Cz in Fig. 7. Therefore, acoustic waves no longer give rise to ill-posedness of PSE for supersonic flows. C,, C,, C, and C, are in the upper half of the t-plane and do not cause ellipticity. Ill-posedness can only arise due to upstream-associated poles in the lower half of the c-plane. We will denote the set of all such poles by S,. We are not aware of any theory concerning the size of S, for compressible flows. If it is infinite and unbounded from below in the t-plane, the inverse Laplace transform cannot be carried out and the problem at hand is still ill-posed in the conventional sense. If S, is either finite (and hence bounded) or infinite, but bounded from below, a Laplace inversion contour can be found below all of these poles, and the requirements for well-posedness are thus satisfied, i.e. the solution of such a system will continuously depend on the initial conditions. However, the upstream exponentially decaying solutions will still be interpreted by the differential operator as downstream exponentially growing ones. Therefore, the parasite solutions will cause contamination. It is in this sense that the problem is still ill-posed. In semi-discretized PSE, the Z-transform can again be applied. As discussed before, the mapping between t-plane and z-plane is given by equation (44) for one-step backward difference. The poles far away from D, pose no problem. A pole such as B, shown in Fig. 7 causes contamination of the solution if it is allowed to enter C,. The minimum step-size required is hence given by equation (47).

4. VERIFICATION

OF THE

THEORY

Although our aim is to study the physical origins of ellipticity in PSE, not numerical methods, we need to test the validity of our analysis via sample computations. Ellipticity due to acoustic waves for incompressible flows gives rise to a step-size restriction (42). This theory is derived and tested by Li and Malik [l 1, 121 and proves to be very accurate despite the frozen-coefficient assumption used in their analysis. This is due to the fact that the eigenfunction of the continuous mode is distributed in the freestream where the coefficients of PSE are constant. When the PSE solution procedure fails, the problem initiates in the freestream. The theory concerning ellipticity due to the discrete upstream mode is tested here by using equation (47) for a compressible boundary-layer flow. We compute the development of a two-dimensional first mode disturbance [22] of frequency F = 0.4 x 10m4in a Mach 1.6 flat-plate boundary layer (F = Sltvf/Uz, f being the frequency in Hz). The spatial eigenspectra are computed at Re = 500 (A’= 500). The complex wavenumber of mode D, is ~1,= 0.04517 - iO.001175 normalized by scales at Re = 500. The upstream mode closest to D, in the sense of equation (48) has wavenumber IX,,= - 0.07395 - iO.07891, whose position in the complex r-plane is approximately given by 5, = u, - al. Then, according to equation (47), the local minimum step-size

Parabolized

stability

equations

295

is found to be 7.7. Farther downstream at Re = 1000 (X = 2000), the minimum step-size drops to approximately 7 (still normalized with scales at Re = 500) due to the slow change of the boundary layer. The PSE computations show that, when 6x = 7 -=z6~,,,~,,the solution is contaminated and oscillations appear in the growth rate curve shown in Fig. 8. When 6x = 8.5 > 6x,,, is used, the PSE solution is smooth. The wavelength of the oscillations is measured from Fig. 8 to be approximately 51. Since Real(5,) x Real(cr, - cr,) = 0; 11912, the wavelength of the parasite solution predicted by theory is Iz = 2n/]Real(Ql = 52.75. This confirms that the ellipticity is caused by this upstream discrete mode. All discrete modes are associated with the high velocity gradients inside the boundary layer. Therefore, contrary to the conclusions of Haj-Hariri [lo] and Herbert [7], the convective terms in PSE are not irrelevant in causing ellipticity. The maximurn in the eigenfunction amplitude of the upstream mode in this example occurs near the wall and hence, when the PSE solution procedure fails, the problem initiates inside the boundary layer. The eigenfunctions of D, (the first mode) and the upstream mode are shown in Fig. 9 for Re = 500. We note that, for the above example, equation (42) yields 6X,i, z 22 while the step-size in the Mach 1.6 boundary layer is governed by equation (47) which yields 6X,i” x 7.7. Hence, the step-size limitation is relaxed for this supersonic case.

5. CONCLUSIONS The difficulties often encountered in implementing the primitive-variable form of PSE have strong physical roots. For subsonic boundary layers, the upstream-associated acoustic branch cut in the spectra of the PSE operator makes the initial-value problem in space ill-posed. In this case, any small perturbation to the initial condition on this branch cut will grow unbounded, and hence, the solutions do not continuously depend on the initial conditions. Another branch cut that was thought, in the previous studies, to represent a viscosity-associated upstream mode is shown to be a downstream convective mode and not responsible for ill-posedness. For supersonic boundary layers, ill-posedness originates from upstream-associated poles. These discrete modes are associated with the high velocity gradient in the boundary layer, and therefore, the convective terms of PSE are important in the upstream propagation of information. The fact that, in practice, PSE can be used successfully is the consequence of the finite-difference model of the original equations. Since the desired solution for PSE is slowly-varying, relatively large discretization step-sizes can be used without destroying the physics which the PSE intends to represent. In the semi-discretized PSE, we have shown that the upstream modes are manageable and pose no threat to the solution

o.oc’12 ti

Fig. 8. Growth rates of the first mode disturbance in a Mach 1.6 flat-plate boundary layer computed using tix < sxrnln and 6,~ > 6x,,,,. Growth rate and x normalized by length scale at Re = 500. CAF 2613-D

by

296

F.

Li and M. R. Malik

(4

-

downstreammode - upstreammode

- 0.6

I. 0

50

I

0,

IUI

1 100

0.6

,,

I,

I 150

I,,

,

200

Phase ( degree ) Fig. 9. Eigenfunctions ic of the downstream first mode disturbance and the upstream mode responsible for ellipticity (Mach 1.6 boundary layer, Re = 500). (a) amplitude and (b) phase.

procedure once the step-size is sufficiently large for which a condition is derived. A computational example is given for a supersonic boundary layer to verify this condition. Finally, it should be pointed out that the present analysis does not deal with the nonlinear iterative process of finding the wave number, tl. This process itself could be subject to numerical instability and is discussed in Li and Malik [ 121.Furthermore, the effect of boundary-layer growth has been ignored in the present analysis. Acknowledgemenrs-The authors wish to thank their colleagues Drs C.-L. Chang and R.-S. Lin for many helpful discussions. This work was sponsored under NASA Contract NASI-20059.

REFERENCES I. Herbert, Th., Boundary-layertransition-analysis and predictionrevisited.AIAA Paper No. 91-0737, 1991. 2. Bertolotti, F. P., Herbert, Th. and Spalart, P. R., Linear and nonlinear stability of the Blasius boundary layer. J. Fluid. Mech., 1992, 242, 44474.

Parabolized stability equations

297

3. Bertolotti, F. P and Herbert, Th., Analysis of the linear stability of compressible boundary layers using PSE. Theorei. Comput. Fluid Dynamics, 1991, 3, 117-124. 4. Chang, C-L., Malik, M. R., Erlebacher, G. and Hussaini, M. Y., Compressible stability of growing boundary layers using parabolized stability equations. AIAA Paper No. 91-1636, 1991. 5. Rubin, S. G., A review of marching procedures for parabolized Navier-Stokes equations. Proc. Symp. on Numerical & Physical Aspects of Aerodynamic Flows. Spring-Verlag, New York, 1981, pp. 171-186. 6. Malik, M. R. and Chang, C.-L., Boundary layer stability and transition. In Handbook of Fluid Dynamics and Fluid Machinery, ed. J. Schetz and A. Fahs. John Wiley, 1996, pp. 318-354. 7. Herbert, Th., F’arabolized stability equations. AGARD-R-793, 1994, pp. 4-1-4-34. 8. Malik, M. R., Li, F. and Chang, C.-L., Crossflow disturbances in three-dimensional boundary layers: nonlinear development, wave interaction and secondary instability. J. Fluid Mech., 1994, u16, l-36. 9. Kreiss, H.-O. and Lorenz, J., Initial-Boundary Value Problems and the Navier-Stokes Equations. Academic Press, New York, 1989. 10. Haj-Hariri, H., Characteristics analysis of the parabolized stability equations. Studies in Appl. Math., 1994, 92,41-53. II. Li, F. and Malik, M. R., Mathematical nature of parabolized stability equations. Proc. of IUTAM Symp. on Laminar-Turbulent Transition, Sendi, Japan, Springer, Berlin-Heidelberg, 1995, pp.202-212. 12. Li, F. and Malik, M. R., On the nature of the PSE approximation. Theoret. Comput. Fluid Dynamics, 8,253-273, 1996. 13. Bertolotti, F. P., PSE for incompressible flows. In Instability & Transition: Experiment, Theory and Computation, ed. T. C. Corke, G. Erlebacher and M. Y. Hussaini. Oxford University Press, (in preparation), 1996. 14. Salwen, H. and Grosch, C. E., The continuous spectrum of the Orr-Sommerfeld equation. Part 2. Eigenfunction expansions. J. Fluid Mech., 1981, 104, 445465. 15. Ashpis, D. E. and Reshotko, E., The vibrating ribbon problem revisited. J. Fluid Mech., 1990, 213, 531-547. 16. Balakumar, P. and Malik, M. R., Discrete modes and continuous spectra in supersonic boundary layers. J. FIuid Mech., 1992, 239, 631-656. 17. Gustavsson, L. H., Initial-value problem for boundary layer flows. Phys. Fluids Al, 1979, 22, 1602-1605. 18. Carrier, G. F. and Pearson, C. E., Partial Differential Equations, Theory and Technique. Academic Press, New York, 1976. 19. Schmid, P. J., Henningson, D. S., Khorrami, M. R. and Malik, M. R., A study on the sensitivity for hydrodynamic stability operators. Theoret. Computat. Fluid Dynamics, 1993, 4, 227-240. 20. Khorrami, M. ‘R. and Malik, M. R., Efficient computation of spatial eigenvalues for hydrodynamic stability analysis. J. Computat. Phys., 1993, 104, 267-272. 21. Sirovich, L., Aa Introduction to Applied Mathematics. Springer-Verlag, New York, 1988. 22. Mack, L. M., Boundary layer linear stability theory. AGARD-R-709, 1984, pp. 3-1-3-81.