On simulation and analysis of instability and transition in high-speed boundary-layer flows

Computing Systems in Engineering, Vol. 6, No. 6, pp. 563-575, 1995

Copyright 0 1995Elsevier Science Ltd Printed in Great Britain. All rights reserved 0956-0521/95-$9.50 + 0.00

Pergamon

ON SIMULATION AND ANALYSIS OF INSTABILITY AND TRANSITION IN HIGH-SPEED BOUNDARY-LAYER FLOWS? C. DAVID F’RUETT~and THOMAS A. ZANG$ $ Program in Applied Science, The College of William and Mary, Williamsburg, VA 23187, U.S.A. 4 Multidisciplinary Design Optimization Branch, NASA Langley Research Center, Hampton, VA 23681-0001, U.S.A. Abstract-The simulation of instabilities and laminar-trubulent transition in high-speed boundary-layer flows represents one of the major computational challenges of the decade. By taking advantage of recent advances in computational science and instability theory for compressible flows, we have formulated an approach to this problem that combines parabolized stability equation (PSE) methodology with spatial direct numerical simulation (DNS). The relatively inexpensive PSE method is used to explore the parameter space, to compute the early (weakly and moderately nonlinear) stages of laminar breakdown, and to provide inflow conditions for the spatial DNS, which is then used to compute the highly nonlinear laminar-breakdown stage. The approach is made feasible by an accurate and efficient DNS algorithm, which has been implemented in parallel on a CRAY C90 supercomputer. The design of the DNS algorithm is discussed in detail, with emphasis on factors that affect both accuracy and efficiency. The method is applied to the investigation of the laminar breakdown of the boundary layer on an axisymmetric sharp cone in Mach 8 flow. Techniques for analysis of the resulting data are also addressed, including novel computational flow imaging procedures.

1. INTRODUCTION By definition, direct numerical simulation (DNS) is the numerical solution of the unsteady Navier-Stokes equations without approximations and without recourse to turbulence or subgrid-scale models. In concept, the fluid motions are resolved down to the Kolmogorov length scale, at which eddies succumb to viscous dissipation. Scaling arguments’ show that the number of grid points for fully resolved DNS scales as Re914, where Re is the Reynolds number. Moreover, the total computational workload scales as Re3. Consequently, the computational requirements for well-resolved DNS are staggering. With the exception of the few highly specialized simulations of turbulence in a periodic box, most direct simulations cannot be considered fully resolved when subjected to rigorous scrutiny. In practice, numerical dissipation in some form usually dampens eddies at scales significantly larger than that of the Kolmogorov scale. Moreover, some level of approximation inevitably occurs; for example, through the imposition of boundary conditions or through empirical viscositytemperature relations. Nevertheless, we retain the terminology “DNS” in a loose sense when applied to unsteady Navier-Stokes simulations that are moderately well resolved. For reasons of computational efficiency, most DNSs have been temporal. In temporal DNS, one TResearch begun under NASA contract NASl-19831 while the first author was employed by Analytical Services & Materials, Inc. and completed under NASl-19656. CSE 6160

examines a window within the flow field that travels downstream at some characteristic velocity, usually that of the phase velocity of a traveling disturbance of interest. Within the window, the flow is assumed to be both locally parallel and periodic in the streamwise spatial dimension. The box length is typically only one or two wavelengths of the fundamental of the disturbance, which contributes to the computational efficiency of the approach. The temporal evolution of the flow within the box and its true spatial evolution are related approximately by a Gaster-like transformation.2 In contrast, in spatial DNS, a flowthrough domain with inflow and outflow boundary conditions is considered. The spatial domain is sufficiently long to incorporate the entire phenomenon of interest (i.e., wave propagation, laminar breakdown, or the complete laminar-turbulent transition process). In general the spatial model is considered to be more physically correct for simulating the natural and forced transition phenomena that occur in wind tunnels and on flight vehicles. Indeed, at least two situations exist in which the temporal ansatz is misleading. First, if the wave propagation phenomena are highly dispersive (disturbances of different frequencies travel at different velocities), then spatial theory is necessary. Second, for compressible wall-bounded flows, the growth of the boundary layer is difficult to account for in a temporal context (although Guo et al.’ have made considerable progress in this area). Because the change of length scale associated with a growing boundary layer can lead to the tuning or detuning of preferred modes of instability, such effects should be incorporated in simulations of transition. Here, our principal 563

564

C. David Pruett and Thomas A. Zang

concern is the physics of transition in compressible flows; consequently, we restrict attention to spatial DNS. Although spatial DNS for incompressible flows is fairly routine at the present time, DNS for compressible flows has received attention only relatively recently. Among the spatial simulations of compressible flows are those of Grinstein et ~f.,~ Lele,’ Thumm,6 Maestrello et al.,’ Eissler and Bestek,8 Pruett et uI.,~ Pruett and Chang,” Kennedy,” Guo a et al.” and Rai et al.‘” Of these simulations, few consider only two-dimensional flows or threedimensional flows with a few spanwise modes. Thus, only a relative handful of spatial simulations of transition and/or turbulence exist for compressible flows. Because of physical and computational considerations, DNS for compressible flows is considerably more daunting than for incompressible flows; this complexity partially explains the small number of practitioners of DNS for the compressible case. Our specific interest is the particularly challenging problem of simulating instability and transition in a spatially evolving high-speed boundary-layer flow, a problem that has been intractable until quite recently. Fortunately, several recent advancements have lead to a new and promising approach to the simulation of transitional high-speed flows. First, the parabolized stability equation (PSE) method has emerged as a new computational and analytical tool that treats both mean-flow nonparallelism and moderately nonlinear wave interactions.‘4m’7 Second, high-order methods of spatial discretization have been adapted to Third, a class of high-order lowspatial DNS. 59~‘n~19 storage methods of explicit time advancement has emerged.20,2’ Fourth, the recent class of mainframe supercomputers (e.g., the CRAY C90) is (nearly) sufficient in speed and memory to accommodate the problem.

_--

Advances in computational science come about by improvements in computer technology, by enhancements in algorithms, and by the judicious matching of the problem, the algorithm, and the machine. Our work falls primarily into the thrid of these categories. We have exploited the recent theoretical, algorithmic, and computational advances mentioned previously to develop a novel approach to the simulation of transition in high-speed flows and to design an accurate and efficient DNS algorithm that is well suited to the problem. Our approach combines the PSE and DNS methodologies in a synergistic manner. The relatively inexpensive PSE method is used to explore the parameter space, to compute the early (weakly and moderately nonlinear) stages of laminar breakdown, and to provide inflow conditions for the spatial DNS, which is then used to compute the highly nonlinear laminar-breakdown stage. The approach is made feasible by an accurate and efficient DNS algorithm that exploits low-storage Runge-KuttaZo methods for time advancement and that incorporates spectral and high-order compact-difference techniques for spatial discretizations. The algorithm is optimized to run in parallel on supercomputers; currently, it executes in excess of 2.5 Gflops on eight processors of a CRAY c90. In the next section, we describe the governing equations and the problem formulation. In Section 3, we discuss a number of algorithmic details, including the trade-off between computational effort and storage, temporal and spatial discretizations, estimation of the numerical stability limit, and the need for filtering and efficient implementation of the filter. In Section 4, we present some results of a simulation of the laminar breakdown of a hypersonic boundarylayer flow along a cone. Novel flow-visualization procedures are also discussed. In the last section, we briefly present some conclusions drawn from this work.

-_(--I-

Fig. 1. Body-fitted

~(x,W

coordinate

-

-

system.

-

-

-

-

Instability 2. PROBLEM

and transition

0

FORMULATION

We consider a three-dimensional boundary-layer flow on an axisymmetric body (Fig. 1) in the bodyfitted coordinate system x = [x, 0, z]‘, where x is the arc length along the body, 0 is the aximuthal angle, and z is the coordinate normal to the body. Associated with x is the fundamental metric tensor g,,, which has the nonzero components g,, = s’;

gz2 = r2;

g,, = 1

-$h2 - (DO,+)& sin $ +pF,+--3 r

H=

s 8x ’

y a0

of,

au 8Z

’

1 a(rsu) 1 i?(m) DIU Dfu =--. I =-----. rs 32 rs 2.x ’

(2)

In the coordinate system of Fig. 1 and the notation of Eq. (2). the dimensionless, compressible NavierStokes equations (CNSE) assume the following conservative form:

aQ

-_=-V

(3)

at

Q = [P>PU, PO, P”, &I’ V=D;E+D;F+DfG+H

PU

r

0

(1)

=_-

cos 4

cos fp -~4+(D’$$)E, r

(6)

K

h,=

where r(x, z) = R + z cos 4, s(x, z) = 1 - z(d$/dx), and R(x) and 4(x) are the body radius and the angle of the surface tangent, respectively (Fig. 1). For notational and computational convenience, we define the following partial differential operators, which incorporate the metric quantities r and s: 1 0’24 I au Dou Y =__. @ju =--.

565

in high-speed boundary-layer flows

(y - 1)M’Re O”,’

h,=

h, =

K (y - I)M’Re

Dir

(y - I;M’Re

of T

where u, u, and w are velocities in the e,. e,), and e, directions, respectively, and E, = p/(y - 1) + (p/2) (u’ + u2 + w’), p, p, T, p, and K are the total energy, density, pressure, temperature, viscosity, and thermal conductivity, respectively. Moreover,

(8) is the stress tensor, 6,, is the Kronecker delta, d = Dtu + Div + Dlw is the divergence of the velocity rate-of(dilatation), and gi, is the symmetric deformation tensor, whose components are

(4)

o,, = Do,u - wDo,qb

(5)

0 22=D;v

cr33 =

+-

u sin f#~ w cos 4 +p r r

Dfu

PUU-511 -P E= (E,

+pju

-

PUU-

712

PUW -

713

UT,1 -

0712 -

g13 = 03, = f(D;w

Pt’ pl?u F=

~~~~~~~~~ 72,

pvv -Ts,,-p pvw (E, + p)v -

UT~, -

723 ~72~ -

~97~~ -

h,

W733 -

h,

PM’

G=

PWU-

731

pwv -

732

pww (E, + p)k” -

U73, -

7j3

- p

u732 -

v sin 4

Do,u+D;v-WT,3 -h,

(7)

1

+

r

D;u - uDo,qb)

D;w +D:v

--

v cos ip r

(9)

With regard to notation, all numerical subscripts refer to physical vector or tensor components, except for Eq. (l), in which the components of the fundamental metric tensor are covariant. The governing system is closed by an equation of state, which is assumed to be that for a perfect gas. With regard to the nondimensionalization of the governing equations, all flow variables except pressure have been scaled by reference values, which we denote with the subscript r. Lengths are scaled by the boundary-layer displacement thickness S * at some

566

C. David Pruett and Thomas A. Zang

reference station x0*, pressure is scaled by p,* u,*‘, and time is scaled by 6 */u: (Throughout this paper, dimensional quantities are denoted by a superscript *.) The dimensionless parameters that arise are the Mach number M, the Reynolds number, Re, the Prandtl number Pr, and the ratio of specific heats y. Here, Pr = 0.7 and y = 1.4. Viscosity and thermal conductivity are each assumed to vary in accordance with Sutherland’s law and are related by K = p/Pr. With proper interpretation, Eqs (l)-(9) are valid for flows over either two-dimensional or axisymmetric bodies. The two-dimensional case is recovered as r - 1 in Eq. (2) and I/r -0 in H [Eq. (6)] and in oli [Eq. (9)]. A nonconservative but computationally useful alternative form of the energy equation is Q,=p;

E,=pu;

F,=pv;

H5 = @ - (y - 1)pd;

G,=pw

@ = Q,Z~

(IO)

where repeated indices imply summations and the subscript 5 refers to the last element of the state and flux vectors. In subsequent discussions, we will refer to Eq. (10) as the “pressure equation” as distinguished from the “energy equation.” Mathematically, DNS is the numerical solution of an initial boundary value problem (IBVP). Consider the computational domain x,” < xi < x0~, O<&

n

0 < rk < z,,,

(i =o,.

(j =0 ,

‘.? N,)

,N,)

(k = 0, . . , N,)

(11)

where n is an integer that specifies the fundamental period in the periodic azimuthal direction. For spatial DNS, inflow and outflow boundary conditions are required along boundaries xi,, and x,,,, respectively. For the Navier-Stokes equations, all five flow variables should be specified at an inflow boundary.” For simulations of instability and transition, the inflow condition is typically comprised of the superposition of a time-dependent (base) state Qo, which is obtained from a spectrally accurate boundary-layer fluctuation Q,, which is code,23 and a time-periodic derived from PSE methodology.‘6 As a theory of fluid instability, PSE treats both nonparallel base states and nonlinear wave interactions. Consequently, Q, can contain a multitude of harmonics of fundamental two- or three-dimensional instability waves. At the outflow boundary, we adapt the buffer-domain concept, 24whereby both the base state and the governing equations are modified to ensure that all waves propagate out of the domain. At the far-field boundboundary ary z = z,,,, we adapt the nonreflecting condition of Thompson,25 which is based on inviscid characteristic theory. At the wall, no-slip conditions are usually imposed on all velocity components. The

thermal boundary condition isothermal or adiabatic.

at the wall can be either

3. ALGORITHM DESIGN

We now consider an algorithm that is particularly appropriate for simulating laminar-turbulent transition at relatively large Mach numbers. For such a problem, temporal and spatial accuracy are overriding considerations. Not only must the method accurately predict the temporal-spatial evolution of primary instability waves, it must further resolve all significant harmonics, which have short wavelengths and high frequencies relative to their fundamental waves. Two primary considerations motivate the use of fully explicit temporal discretizations. First, because /J and K vary considerably with the large temperature fluctuations associated with high-speed wall-bounded flows, the diffusion terms of the CNSE are moderately nonlinear. Thus, the efficient semi-implicit techniques developed for the incompressible Navier-Stokes equations are not immediately applicable. For the simulation of boundary-layer phenomena in high-speed compressible flows, the advection and diffusion constraints on the time step may be of the same order of magnitude. Hence, no computational advantage is gained by the removal of the diffusion constraint by semi-implicit methods. Second, although fully implicit methods permit a larger time step, the gain is offset by additional storage requirements and by the necessity of iteration within each time step. Furthermore, the time step is more constrained by considerations of the accuracy needed for the high-frequency harmonics than by stability. Thus, we have opted for fully explicit temporal discretization using low-storage technology. Here, we describe a simple but highly accurate algorithm that is built on fully explicit temporal discretization, nondissipative high-order methods of spatial discretization, and modest filtering. Spatial discretizations In simulations of instability and transition, the growth rates of instability waves must be computed accurately; for this purpose, we prefer nondissipative methods. Consequently, we exploit spectral and centered compact-difference methods. Here, we describe a collocation method based on these discretizations Let the computational domain [Eq. (1 I)] be discretized into N,, N,, and N, points in the x, 0, and z coordinate directions, respectively. Among the discrete derivative operators we have considered are the following: 1. Streamwise a. compact-difference b. compact-difference c. compact-difference et al.26)

3-4-3 scheme (Pade) 3,4-6-4,3 scheme (Lele’) 5,5-6-&S scheme (Carpenter

Instability

and transition

in high-speed

2. Spanwise

(azimuthal) a. Fourier spectral collocation with asymmetry allowed (Canuto et aL2’) b. spectral collocation with sine and cosine series (symmetry enforced) 3. Wall normal a. Chebyshev spectral collocation (Canuto et ~1.~‘) b. compact-difference 3-4-3 scheme (Pade) c. compact-difference 3,4-6-4,3 scheme (Lele’) d. compact-difference 5,5-6-55 scheme (Carpenter et ~21.~~)

In regard to the compact-difference schemes, the notation 3,4-6-4,3, which refers to the method of Lele,’ for example, indicates that the differentiation operator is formally sixth-order accurate at the interior of the stencil, with boundary closures that lead to third-order accuracy at boundary points and fourth-order accuracy at points adjacent to the boundary. Readers are referred to the cited papers for details of the specific methods. Unlike conventional finite-difference operators, compact-difference and spectral operators are global. To allow for this multiplicity of differentiation options in a transparent manner, spatial derivative evaluations are accomplished by subroutine calls with the following syntax: CALL PARTIAL (isym,m,u,du,iopt). Given the flow-field quantity “u”, subroutine PARTIAL returns in “du” the partial derivative of “u” with respect to x,. The “isym” option, which is operative whenever m = 2, enforces azimuthal symmetry if isym = 1, which reduces the computational effort by a factor of 2 with the use of fast sine and cosine transforms. Otherwise, complex Fourier transforms (FFTs) are exploited. The option “iopt” controls whether or not metric terms are incorporated into the operator. For example, for u = u and m = 1, D~ZJ is returned if iopt = 1; otherwise, D(lu is returned. By this convention, the routine that evaluates the right-hand side of Eq. (3) is considerably simplified, and the same code readily treats either axisymmetric or two-dimensional geometries. Also passed to PARTIAL (via common block) are the parameters “lmethx,” “lmethy”, and “lmethz”, which select the differentiation operators in the streamwise, azimuthal, and wall-normal dimensions, respectively, in accordance with the list above. The results presented here were obtained with the 5,5-665,5 and 3,4-6-4,3 compact-difference schemes in the streamwise and wall-normal directions, respectively, and the symmetric sine and cosine transform method in the periodic (azimuthal) direction. For reasons we cannot yet explain, we have been unable to maintain numerical stability with the method of Carpenter et a1.26 in both the streamwise and wall-normal directions. Although the Chebyshev spectral collocation method in the wall-normal direction is useful for certain diagnostic purposes, it imposes an excessively small constraint on the time step, which renders the method impractical for large-time simulations.

boundary-layer

561

flows

Table 1. Minimum number of derivative evaluations in assembly of V in Eq. (3) No.

Calls to ‘PARTIAL Rate-of-deformation tensor Temperature gradient VT Flux evaluations (3 x 5) Total derivative

In

direct

9 3 15

umn

21

calls

simulation,

there

is a design

trade-off

CPU requirements. Most of the computational effort is devoted to computing the partial derivatives of flow-field quantities with the methods described above. For example, with fast transform techniques, the evaluation of derivatives with respect to 0 at all points in the flow field requires Compact-difference 0 (N, N, N, log, No) operations. derivatives, which involve the solution of multiple tridiagonal systems of equations, require O(N, Ni N,) operations. Given that the time-advancement scheme is fully explicit, the computational effort per time step must be minimized. The minimum number of partialderivative calls needed to evaluate the compressible Navier-Stokes equations is 27, as shown in Table 1. Because we favor minimal work over minimal storage, we construct the algorithm so that only the minimum number of calls to PARTIAL is required. between

memory

and

Temporal discretization The spatial discretization and evaluation of V described above reduces Eq. (3) to a system of ordinary differential equations of order N, x N, x N,. We now seek a time-advancement scheme for which storage is a minimum, given the minimal-work constraint above. Table 2 summarizes the minimum number of threedimensional arrays required when a third-order low-storage RungeKutta (RK3) scheme2’ is used to advance the time. Table 2 shows clearly why low-storage RK methods are so advantageous. A standard fourth-order Runge-Kutta (RK4) method, for example, would require 20 arrays for evaluation of the right-hand side, which effectively doubles the total storage requirement (and renders the present computation impractical on existing supercomputers). Note that use of the pressure rather than the energy equation saves two temporary storage arrays because evaluation of the dissipation function @ requires less temporary storage than does evaluation of energy fluxes. Although the pressure equation is nonconservative, we prefer it for this reason. Table 2. Minimum storage requirements for low-storage Runge-Kutta Arrays

of dimension

time advancement

O(N, x N, x Nz)

Primitive variables @, U, a, W,p, p) Low-storage Runge-Kutta right-hand side Temporary work space of pressure (energy) Total

of pressure

(energy)

form

No.

form

6 5 5(7) 16(18)

568

C. David

Pruett

and Thomas

If fully explicit schemes are to be useful in practice, the numerical constraint on the time step must be estimated accurately. Consequently, we adjust the time step dynamically so that we continually maintain nearly maximal time increments. Specifically, we require

JliGZ.A,

A. Zang RK3 2.5

stability

region

I

I

-2

-1


lim,)

(12)

where 0 < sf < 1 is the stability limit safety factor (typically sf = 0.98), and where lim,, lim,, and lim, are limiting values obtained by independent consideration of the stability of model linear advection, viscous diffusion, and thermal diffusion equations, respectively. For example, to evaluate the advection limit, we consider the linear transport equation q,+(lul+c)q,+(l~I+c)q,,+(l~~‘I+C)ql=O

lim, < 1

(14)

lim, =

(l~I+c),k+(lt’l+~),k+(/~~l+c),,k CFL;Az,

CFL,. Ay,

1

(15)

The values CFL,, CFL,., and CFL,X are given in Table 3; these values were obtained by numerical experimentation or by Von Neumann stability analysis. Similarly, lim, and lim, are defined by considering model diffusion equations of the form

q&V2q;

VE

f,” {P P

I

(16)

with the help of Table 4.

Table 3. CouranttFriedrichs-Levy (CFL) limit 1c(A,t /Ax)) < CFL for scalar linear advection equation 4, + “4, = 0 Stability limit CFL

Spatial operator Compact 4th-order Compact 6th-order Fourier spectral Chebyshev spectral

Table 4. Viscous

1.0

Method of derivation Empirical Analytical Analytical Empirical

stability limit Iv(At/Ax’)1 < VI for scalar diffusion equation 4, = ~4,~ Stability

Spatial operator Compact 4th-order Compact 6th-order Fourier spectral Chebyshev spectral

limit

VI

2.5113 2.5114 2.51/r? 2(2.51/n*)

0

diffusion Fig. 2. Stability domain of third-order Runge-Kutta with embedded ellipse. Stability limits of test problem by “x”.

method marked

(13)

where velocities u, v, and w and the speed of sound c are obtained from the true problem. The stability limit of Eq. (13) can be shown to be At

-3

Method of derivation Empirical Empirical Analytical Empirical

The domain of absolute stability for RK3 methods (Fig. 2) can be found in many references, including Canuto et aL2’ The factors of 2.51 and & that arise in Tables 3 and 4, respectively, pertain to the intercepts of the RK3 stability boundary with the real and imaginary axes, respectively. The stability of a pure diffusion problem is limited by the value -2.51 of the real-axis intercept, whereas the stability of a pure advection problem is limited by the values k ,,/? of the imaginary-axis intercepts. We conjecture that solutions of problems of mixed advection-diffusion type remain stable provided, for example, that the point (-2.51 lim,, fi. lim,) lies with the RK3 domain of absolute stability. Equation (12) is a convenient approximation to this criterion based on the fact that a half-ellipse through the real and imaginary intercepts lies entirely within this region. The point identified by the “xl’ in Fig. 2 identifies limiting values that are typical of the actual test case of Section 4. The numerical stability of semidiscrete and fully discrete compact-difference schemes is subtle and is currently an area of active research.26 In general, what is true for scalar equations is not necessarily true for systems. Our method of estimating stability limits is somewhat heuristic but has proved to be a useful guideline in practical applications for a wide variety of numerical methods and over a wide range of parameter values. For the forced instability problems that we have considered, 700-1000 time steps per period of the fundamental disturbance suffice to maintain numerical stability. Recently, Carpenter and Kennedy” have derived a family of fourth-order, five-stage low-storage RK schemes with a considerably extended region of absolute stability. The increased effort of two additional stages relative to the RK3 scheme is more than compensated for by a significant increase in time-step size. These new schemes promise a 1 l-16% increase in computational efficiency, with higher accuracy

569

Instability and transition in high-speed boundary-layer flows 1.2

I _--__

I

I

’

I

’

I

(

Carpenter

1 .o

0.8 0.6

0.8 zk 8 ,-? 0.6

0.4 0.2

+ -

0.0

0.4

aAx

0.2 Fig. 4. Transfer

0.0

0.0

0.2

0.6

0.4

0.8

low-pass

compact-

1 .O

aAx n Fig.

function for sixth-order difference filter.

L

3. Effective wave number versus wave number selected compact-difference schemes.

for

as a bonus. These schemes have been successfully implemented for temporal DNS of compressible flows in unpublished work by Erlebacher (personal communication). Filtering

The eigenvalues of first-order central-difference operators (denoted by D) are known to vanish for periodic spatial modes exp(lax) (where I = 0) such that aAx = x, as shown in Fig. 3. For reasons of computational efficiency, we prefer to evaluate the viscous terms of the form D( pDu)

= DpDu

+ pD2u

(17)

by using the left hand side of Eq. (17). An unintended effect is that no viscous damping is applied to the “sawtooth” mode, the result of which for nonlinear problems is nonphysical two-point oscillations in the solution. Several possible fixes exist, including upwinding, staggered grids, evaluation of viscous terms with the right-hand side of Eq. (17) (where D2 is a second-derivative operator), and filtering. From a practical point of view, filtering is the simplest approach and perhaps the most versatile. In our algorithm, we apply a low-pass sixth-order compactdifference filters to the solution every few time steps. The filter has two free parameters, which we use to adjust the amount of dissipation. For our purposes, we desire minimal dissipation. Typically, we use parameter values for which the coefficient of the sixthorder truncation error term is $(A.x)~; the corresponding transfer function is shown in Fig. 4. We find that filtering is necessary only in the wall-normal direction and only approximately every third full time step. When amortized over the total computational

effort, the additional effort of filtering is very small (less than 2% of the total effort). The elimination of numerical high-frequency oscillations also virtually eliminates nonphysical numerical reflections from the outflow boundary. (See Ref. 9.) Vectorization

and parallelization

Finally, we add a few comments in regard to parallel implementation of the algorithm. The algorithm as originally structured was nearly the perfect candidate for autotasking on CRAY supcrcomputers. For example, in most subroutines, inner loops range over (x, f3) surfaces; thus, vectors are typically of length (N, + 1) x (N,) + 1). Parallel outer loops typically range over index 0 < k < N,. For certain subroutines (e.g., tridiagonal solvers and Fourier transforms), recursions in the outer loop inhibit parallelization, in which case we “strip mine” to take advantage of both parallelization and vectorization within the inner computational loops. The algorithm currently runs at 420 Mflops on a single processor of the CRAY C90 and at 2.5 Gflops on eight processors, for a speedup of 5.9. Little effort was spent on optimization, and almost certainly these already respectable performance factors could be improved significantly. 4. FLOW ANALYSIS AND

VlSUALlZATlON

Spatial DNS is the computational analog of forced instability and transition experiments, and the practitioner of spatial DNS often refers to him/herself experimentalist”. Consequently, as a “numerical many of the most useful analysis and visualization techniques for DNS are motivated by their counterparts in physical experiments (e.g., Fourier timeseries analysis). An advantage of simulation over experiment is that the former provides a level of flow-field detail unattainable by any other technique. Conversely, however, the amount of data that is generated by spatial DNS is far greater than current storage allotments allow to be saved or processed. Hence, choices and compromises are necessary in

570

C.

David Pruett and Thomas A. Zang

postprocessing DNS data. In this section, we describe our approach to flow analysis and flow visualization for a specific numerical experiment.

The relationship between Re, and the conventional Reynolds number Re, (which is based on x*) is

Test case

Flow analysis

The test case is the approximate computational analog to the experiment of Stetson et al.,28 which investigated instability and transition on an axisymmetric sharp cone in Mach 8 flow. Specifically, the cone has a half-angle of 7”, which results in a shock angle of 10.4” and a post-shock Mach number of 6.8. The free-stream unit Reynolds number is 1.0 million per foot, which yields a post-shock unit Reynolds number of 1.433 million per foot. For the purposes of the simulation, flow conditions are considered to be constant along the edge of the boundary layer. For high-speed wall-bounded flows, multiple instability modes coexist (see the pioneering work of Mackz9), of which the first and second modes are most important. In the physical experiment, the dominant instability waves are of the second-mode type and lie in the frequency range of 70-150 kHz. Considerable information is available in regard to a particular disturbance that corresponds to a frequency of 102 kHz. From this information, a viable transition scenario was reconstructed by Chang and Malik,” who performed parameter studies with nonlinear PSE methodology. The PSE method is a new stability analysis tool that treats both mean-flow nonparallelism and moderately nonlinear wave interactions. These studies suggested oblique second-mode breakdown as a likely path to transition ” for the geometry and flow conditions of interest. In this scenario, the nonlinear interaction of a pair of equal and opposite secondmode instability waves generates strong streamwise vorticity, which eventually leads to laminar breakdown. (A more thorough discussion of oblique second-mode breakdown can be found in Refs 10 and 17.) After a viable transition mechanism was identified, a synergistic approach that combined PSE and DNS methodologies was used for the simulation. Specifically, the PSE method was used to compute the weakly and moderately nonlinear stages of evolution of the imposed disturbances and thereby to provide a harmonically rich inflow condition to the DNS. The strongly nonlinear and laminar breakdown stages were subsequently computed by well-resolved spatial DNS. Details in regard to the parameter values, the imposed disturbances, and the computational domain can be found in Pruett and Chang.” Here, we confine ourselves to analysis and visualization of the DNS results. For consistency with linear stability theory and experimental results, in the presentation of our numerical results to follow streamwise locations are identified by the Reynolds number Re,, which is based on the boundary-layer length scale

One of the most useful analytical tools in physical experiments of instability and transition is Fourier time-series analysis. In numerical experiments, we also rely heavily on time-series analysis; however, practical differences exist between analytical techniques for physical and numerical experiments. In wind-tunnel stability experiments, time-series data are collected by hot-wire anemometry. Practical considerations limit the placement of the hot wire(s) to relatively few locations, at which the recorded time traces are usually quite long. In principle, in DNS one can think of each grid point as having its own hot wire and associated time trace. In practice, however, for a wellresolved simulation, it is feasible to record only a small fraction of the available data. For example, for the present simulation N, = 4320, N, = 48, and N, = 128; thus, more than 27 million grid points were used. Consequently, in spatial DNS we record time traces for two somewhat different purposes. One use of time traces is to assess transient effects. In a forced instability experiment, the initial condition is that of an unperturbed laminar flow into which instability waves propagate. Normally, we are most interested in the fully developed unsteady flow; transient effects are of secondary interest. To assess transient effects, we rely on a small number of time signals obtained at selected points in X. The trace of one such signal (at Re, = 2606) is shown in Fig. 5. From the signal, we can deduce that the fluid at the point is unperturbed from 0 < t < 90, experiences transients associated with the leading edge of the propagating wave front from 90 < t < 120, and is aperiodic (but essentially fully developed) for 120 < t. Aperiodicity at the later times is more clearly seen in the artificial phase portrait of the signal, which we have omitted here. After the flow has settled into its fully developed state, we record the time history of the flow for a few periods at a relatively dense spatial resolution. For

(18)

Re, = 6.

1 .o

z X

0.8 0.6 0.4

80

100

120

140

Fig. 5. Time trace of total density

160

at Re, = 2606.

180

Instability and transition in high-speed boundary-layer flows

5.0

skin-friction x1o-4

shape

coefficient 35

I

*

I

n

571

factor I

I

’

I

I

n

33 31 29 H

27 25 23 21

3.5 2.2

2.3

2.4

2.5

2.6

2.7

2.8 X10”

2.2

2.3

2.4

2.5

2.6

2.7

2.8 X103

Re,

9%

Fig. 6. Mean skin-friction coefficient and mean shape factor.

data collection, we typically include all grid points in azimuth but thin resolution in the wall-normal and streamwise dimensions by factors of 2 and 24-48, respectively. Whereas the flow changes rapidly with x in physical space, the streamwise evolution in Fourier space of various harmonics is relatively gradual; thus, we perform most of the necessary data compression by thinning in the streamwise direction. Typically, we sample at a rate of 64 samples per period of the fundamental disturbance; however, any sample rate can be specified. Logic is included in the timeadvancement scheme to anticipate when a time sample is due and to automatically adjust the RK3 timestep accordingly, so that time samples occur at precise (to machine precision) intervals. As they are generated, time samples are archived for postprocessing. Among the quantities that we currently extract during postprocessing are mean quantities (such as skin friction and shape factor), harmonic energies, Reynolds stresses, turbulent kinetic energy, and dissipation rate. Figure 6 serves as a road map to the laminar-breakdown process by showing the evolution with respect to x of the mean skin-friction coefficient cr and the mean shape factor H. In a laminar boundary-layer flow, the shape factor, which is the ratio of the displacement thickness to the momentum thicknesses, remains virtually constant. However, as the boundary layer transitions, the momentum thickness grows more rapidly than the displacement thickness, and the shape factor falls precipitously. A transitional state is also indicated by changes in the mean u velocity, which produces a dramatic rise in skin friction. Figure 6 suggests that the flow can be considered fully developed for Re, < 2600; however, downstream of this point remnants of the inital laminar state can be found (e.g., the final upturn in H). Figures 6-l 1 were obtained from the time samples of the four-period interval 56 < t, f 60 (where tp is the time scaled by the fundamental period). Data were averaged both in time and in azimuth; subsequently, we will denote averaged quantities by an overline.

The circles shown in Fig. 6 correspond to particular streamwise stations at which details of the flow will be examined. Second-mode disturbances are dominated by temperature and density fluctuations, which attain maxima near the critical layer (at about one displacement thickness from the wall). This fact is reflected in the root-mean-square (rms) density fluctuation shown in Fig. 7. [Throughout this work, primes denote fluctuations about the mean (i.e., p’ = p - p).] As laminar breakdown occurs, the region of large fluctuations broadens, and the peak shifts toward the wall. Of particular interest to transition modelers are the (density-weighted) Reynolds stresses P uj UJ Figure 8 presents the evolution of the principal components of the Reynolds stress tensor. We observe that the Reynolds stresses are dominated by the streamwise component and that all components of the Reynolds stresses exhibit peaks near the critical layer, which is in keeping with the character of second-mode instability waves. The dominance of the streamwise

rms n

density

fluctutation

I

I

I

0.5

7

Prms Fig. 7. Root-mean-square

density

fluctuation.

572

C. David

Pruett

and Thomas

A. Zang

turbulent

kinetic

I

2.0

r

energy

,

8

K

,

v

1.5 Z 6*

-. 1.0

‘...

---__

___----

0.5

0.0

I

4 x1o-3

3

2

1

0

I

I

0.5

0.0

----____ --. _--

_____------

I

I

i

I

1.5

1.0

2.0 x1o-3

PK

pvv

Fig. 9. Density-weighted turbulent kinetic energy PK. 2.0

&I'

v

'

n

'

I-

1.5

dissipation

Z ;sTi

rate

1.0 2.0 0.5 1.5 0.0

Z

0.0

1.0

0.5

1.5

2.0 x1o-4

6’

0.5

Pw’w’ I

1

I

1

I

’

I

1

0.2

0.0

i :

c__

I

I

I

0.6

0.4

I

6

0.8

1.0 x1o-5

P

-------------_________________ -~-~-.-.-._._____._._____

Z a

I

0.0

2.0

1.5

.;, --..-----____ --_, ______.__.... ....-.. ~_._._._._.-.~_______--____.-.-- ____--

1.0

Fig. 10. Density-weighted dissipation rate Pt.

.-. __-------2.

1.0

: .’I’ 0.5

;: ;___...~ I:

0.0

0.2

0.4

0.6

1.0 x1o-4

0.8

Fig. 8. Principal components of density-weighted Reynolds stress: (a) p u’u’, (b) p v’u’, (c)p w’w’.

component of Reynolds stress is reflected in the turbulent kinetic energy K (Fig. 9) which is defined as one-half the trace of the Reynolds-stress tensor pro[i.e., PK = P/~(u:u~)]. Note that corresponding files in Figs 8(a) and 9 are nearly identical in shape. Of particular importance in two-equation turbulence models is the rate at which turbulent kinetic energy is dissipated. The dissipation rate is denoted by t and is defined by pc = @’ = 0;7:,

We present

the evolution

of t in Fig.

(19) 10.

10 -1 $ 3 -c E” 1o-z m p

1o-3

_~~-~~~~~-~~IS‘7r-,~~~~~~.” .3__ . j’4 ,/ -.__ _ ,I ‘ii,;: ,,..“/--__ -L__ .&~._~ip$.-. .:.,’ d$5’ ! i!p&,, ;r+ _/-.

.. ’ ..’ *.. ,,: I. ,_.-..v .__,;8

_

,,..-Y, ,I *o-4

I ;

1.6

I

I

1.8

,i’ :, : !/ I

2.0

I

I

2.2

I

.._.._ --_

:;

____ ._.__ 1: 02 FJ

_

I

2.4

2.6

2.8

x103

ReL Fig. 11. Evolution of temperature maxima of selected harmonics. (First and second integers in legend identify harmonic

with respect

to frequency respectively.)

w and wavenumber

fi,

573

Instability and transition in high-speed boundary-layer flows Finally, if a Fourier transform is applied to the data both in time and in azimuth, we can decompose the evolution of the transitional flow into its many harmonics, a few of which are shown in Fig. 11. Harmonics are each labeled with a pair of integers (m, n), where m and n identify the harmonic with respect to the fundamental temporal frequency w and the fundamental spanwise wave number /I, respectively. For the test case, transition is triggered initially by a fundamental triad comprised of two mirrorimage oblique waves [(l, 1) and (1, - 1)] and one axisymmetric wave (l,O), each with a dimensional frequency of 102 kHz. All other harmonics arise from nonlinear (sum and difference) interactions. Figure 1I shows the evolution of selected harmonics during the early stages of laminar breakdown. This stage is characterized by saturation of the fundamental waves and rapid growth of higher harmonics. Of particular

Fig. 12. Wall-normal

density

gradient

in selected

importance is the (0,2) harmonic, which dominates all others at Re, > 2400. The counterpart of this harmonic also dominated in the temporal simulations of Pruett and Zang18 and Sandham and Reynolds.” In the former of these papers, a numerical experiment was performed in which the (0,2) harmonic was artificially removed; its suppression inhibited transition. Harmonics of zero frequency appear as structures (vortices) aligned in the streamwise direction, which introduce strong spanwise variations into the flow field. Mounting evidence indicates that such harmonics are essential to laminar breakdown, which is fundamentally a three-dimensional phenomenon. Curiously, the rapid changes in shape factor and skin friction in Fig. 6 correspond spatially to the dominance of the (0,2) harmonic in Fig. 11. Whether or not this warrants further coincidence implies causality investigation.

azimuthal

planes:

(a) just downstream

of Re, = 2257

and (b) just upstream of Re, = 2728. (Flow is left to right. Surface of cone is shown in grey.)

514 Flow

C. David Pruett and Thomas A. Zang visualization

A striking feature of pretransitional high-speed boundary-layer flows is the appearance of so-called ropelike structures in schlieren images. (See the excellent review on this subject by Smith.“) These flow features have not been explained completely; however, they are believed to be related to second-mode instability waves and to be made visible because of the large density fluctuations associated with secondmode disturbances. These structures have a characteristic wavelength of approximately two boundary-layer displacement thicknesses, which corresponds to the expected wavelength of the most unstable secondmode disturbances. To investigate the origin of these structures, we examine the wall-normal density gradient ap/az, the motivation being the schlieren images are sensitive to the first derivative of density in some preferred direction. To accomplish this, we use the Flow Analysis Software Toolkit (FAST) developed at NASA Ames Research Center. FAST incorporates an internal calculator, which we utilize to compute the instantaneous gradient of density and to select the wall-normal component of the gradient. Color contours of the density gradient are then displayed in particular (x, z) planes. The density gradient is remarkably rope-like in appearance in certain planes (Fig. 12). Figure 12 also shows considerable variation in different azimuthal planes, and the integrated effect over all planes is difficult to predict. New developments in the field of computational flow imaging (CFI)32 promise that realistic schlieren imaging will soon be available for numerically generated data (Havener, personal communication). This capability should help to resolve the remaining controversy that surrounds the origin of the ropelike structures.

5. CONCLUSIONS

The computational challenge of simulating laminarturbulent transition in high-speed wall-bounded flows is feasible on existing supercomputers, provided that the algorithm is highly tailored to the physical problem. In this work, we have designed a highly accurate algorithm that is based on low-storage Runge-Kutta methods of time advancement and high-order spectral and compact-difference methods for spatial discretization. The algorithm is efficient and runs in excess of 2.5 Gflops on eight processors of a CRAY C90. The method has been used to examine the laminar breakdown of the boundary layer on an axisymmetric sharp cone in Mach 8 flow. The level of flow-field detail obtained from the well-resolved spatial direct numerical simulation (DNS) is currently unavailable from any other method, experimental or computational. The analysis and visualization efforts with this data have begun to answer fundamental questions in regard to transition physics for high-speed flows. In particular, the DNS has confirmed that oblique second-mode breakdown” is a viable path to transition for the

geometry and flow conditions of the present problem (although it may not be the path preferred by nature). Novel flow-visualization procedures, including numerically generated schlieren-like images, have demonstrated that the ropelike waves seen in schlieren photographs from physical experiments are likely to be manifestations of large-amplitude second-mode instability waves, as suspected. Finally, for the early stages of laminar breakdown, Reynolds stresses, turbulent kinetic energy, and turbulent dissipation rates have been extracted from the numerical data. These quantities should be of considerable benefit to transition modeling efforts for high-speed flows.

REFERENCES 1. S. A. Orszag, “Renormalization group theory,” in ICASE/LaRC Short Course on Turbulent Flow Modeling and Prediction, 14-18 March 1994. 2. M. Gaster, “A note on the relation between temporallyincreasing and spatially-increasing disturbances in hydrodynamic stability,” J. Fluid Mech. 14, 222-224 (1962). 3. Y. Guo, N. A. Adams and L. Kleiser, “Modeling of nonparallel effects in temporal direct numerical simulations of compressible boundary-layer transition,” Theoret. Comput. Fluid Dynamics 1996 (to appear). 4. F. F. Grinstein, E. S. Oran and F. Hussain, “Threedimensional numerical simulation of a compressible, spatially evolving mixing layer,” AIAA Pap. No. 888 0042, 1988. 5. S. K. Lele, “Direct numerical simulation of compressible free shear flows,” AIAA Pav. 89-0374, 1989. 6. A. Thumm, “Numerische Uniersuchungen zum laminarturbulenten Stroemungsumschlag in transsonischen Grenzschictstroemuneen.” Ph.D. thesis. Universitaet Stuttgart, 1991. A. Bayliss and R. Krishnan, “On the 7. L. Maestrello, interaction between first- and second-mode waves in a supersonic boundary layer,” Physics of Fluids A, 3(12), 3014-3020 (1991). 8. W. Eissler and H. Bestek, “Spatial numerical simulations of nonlinear transition phenomena in supersonic boundary layers,” In Transitional and Turbulent Comoressible Flows (Edited bv L. D. Kral and T. A. Zana). ‘FED-Vol. 151, pp. 69-76. ASME, New York, 1993: 9. C. D. Pruett, T. A. Zang, C.-L. Chang and M. H. “Spatial direct numerical simulation of Carpenter, high-speed boundary-layer flows-Part I: algorithmic considerations and validation,” Theoret. Comput. Fluid Dynamics, I, 49-76 (1995). 10. C. D. Pruett and C.-L. Chang, “Spatial direct numerical simulation of high-speed boundary-layer flows-Part II: Transition on a cone in Mach 8 flow,” Theoret. Comp. Fluid Dynamics 7(5), 397-424 (1995). 11. C. A. Kennedy, “The numerical simulation of variabledensity compressible shear layers,” Ph.D. thesis, University of California at San Diego, 1994. 12. Y. Guo, N. A. Adams and L. Kleiser, “Direct numerical simulation of transition in a spatially growing compressible boundary layer using a-new Fourier method,” Workshoo on Direct in Proceedinas 1st ERCOFTAC and Large-E;dy Simulations, Kluwer, Boston, 1994. “Direct 13. M. R. Rai, T. B. Gatski and G. Erlebacher, simulation of spatially evolving compressible turbulent boundary layers,” AIAA Pap. No. 95-0583, 1995. “Linear and nonlinear stability of 14. F. P. Bertolotti, boundary layers with streamwise varying properties,” Ph.D. thesis, Ohio State University, 1990.

Instability

and transition

in high-speed

15. F. Bertolotti and T. Herbert, “Analysis of the linear stability of compressible boundary layers using the PSE,” Theo&. Cornpit. Fluid Dynamics j, 117-124 (1991). 16. C.-L. Charm. M. R. Malik. G. Erlebacher and M. Y. Hussaini, “Cbmpressible stability of growing boundary layers using parabolized stability equations,” AIAA Pap. 91-1636, 1991. 17. C.-L. Chang and M. R. Malik, “Non-parallel stability of compressible boundary layers,” AIAA Pap. 93-2912, 1993. 18. C. D. Pruett and T. A. Zang, “Direct numerical simulation of laminar breakdown in high-speed, axisymmetric boundary layers,” Theoret. Comput. Fluid Dynamics 3, 345-367 (1992). 19. N. A. Adams and L. Kleiser, “Numerical simulation of transition in a compressible flat-plate boundary layer,” in Transitional and Turbulent Comwessible Flows-1993 (Edited by L. D. Kral and T. A. Zang). FED-Vol. 151, pp. 101-110. ASME, New York, 1993. 20. J. H. Williamson. “Low-storaee RunaeKutta schemes.” J. Compur. Phys. i5, 48-56 (1580). 21. M. H. Carpenter and C. A. Kennedy, “Fourth-order SIAM J. Sci. 2n-storage Runge-Kutta schemes,” Comp., 1996 (to appear). 22. T. J. Poinsot and S. K. Lele, “Boundary conditions for direct simulations of compressible viscous flows,” J. Comput. Phys. 101, 1044129 (1992). 23. C. D. Pruett and C. L. Streett, “A spectral collocation method for compressible, non-similar boundary layers,” Int. J. Numer. Meth. Fluids, 13(6), 713-137 (1991).

boundary-layer

flows

515

24. C. L. Streett and M. G. Macaraeg, ‘Spectral multidomain for large-scale fluid dynamic simulations,” Applied Numerical Mathematics, 6, 123-139 (1989/90). 25. K. W. Thompson, “Time dependent boundary conditions for hyperbolic systems,” J. Comput. Phys., 68, 1-24 (1987). 26. M. H. Carpenter, D. Gottlieb and S. Abarbanel, “The stability of numerical boundary treatments for compact high-order finite-difference schemes,” J. Comp. Phys., 108(2), 272-295 (1993). 27. C. Canuto, A. Quateroni, M. Y. Hussaini and T. A. Zang, Spectral Methods in Fluid Dynamics. SpringerVerlag, Berlin, 1988. 28. K. F. Stetson, E. R. Thompson, J. C. Donaldson and L. G. Siler, “Laminar boundary layer stability experiments on a cone at Mach 8. Part- 1: Sharp cone,” AIAA Pan. 83-1761. 1983. 29. L. M. M&k, “Boundary-layer linear stability theory,” in Special Course on Stability and Transition of Laminar Flow, (Edited by R. Michel), AGARD Report No. 709, pp. 3.1, 3.81, 1984. 30. N. D. Sandham and W. C. Reynolds, “Threedimensional simulations of large eddies in the compressible mixing layer,” J. FluidMech. 224, 133-158 (1991). 31. L. G. Smith, “Pulsed-laser schlieren visualization of hypersonic boundary-layer instability waves,” AIAA Pap. 94-2639, 1994. 32. G. Havener, “Computational flow imaging: fundamentals and history,” AIAA Pap. No. 94-2615, 1994.