A11107 3“10Cm0 PUBLICATIONS

NBSIR 83-2686

A Numerical Procedure for An

Inviscid Stability Analysis of An

Axisymmetric Jet

U.S. DEPARTMENT OF COMMERCE

National Bureau of Standards

Center for Chemical Engineering

Chemical Process Metrology Division

Washington, DC 20234

May 1 983

999Z-fP J S' DEPARTMENT OF COMMERCE

Qcn * national bureau of standards

NBSIR 83-2686

A NUMERICAL PROCEDURE FOR AN

INVISCID STABILITY ANALYSIS OF AN

AXISYMMETRIC JET

NATIONAL BUREAU

OF STANDARDS

LIBRARY JUN 2 7 1985

C\o-V cccc. - £_a rc. . c3l v oo

g^-d>w8W

Elizabeth F. Moore

U.S. DEPARTMENT OF COMMERCE

National Bureau of Standards

Center for Chemical Engineering

Chemical Process Metrology Division

Washington, DC 20234

May 1 983

U.S. DEPARTMENT OF COMMERCE, Malcolm Baldrige, Secretary

NATIONAL BUREAU OF STANDARDS. Ernest Ambler. Director

A NUMERICAL PROCEDURE FOR AN INVISCID STABILITY

ANALYSIS OF AN AXISYMMETRIC JET

Elizabeth F. Moore

Chemical Process Metrology Division

Center for Chemical Engineering

National Engineering Laboratory

National Bureau of Standards

Washington, D. C. 20234

Inviscid stability analysis has been applied to the mixing

layer profile of an axisymmetric jet and a coflowing stream. A

collection of computer subprograms has been developed to solve

the resulting eigenvalue problem. The effect of changing the

velocity profile and its parameters can be easily assessed.

Results for Gaussian profiles are included.

Key Words: Axisymmetric jet; eigenvalue problem; mixing layer;

numerical methods for eigenvalue problems; stability

analysis.

Introduction

The mixing region between a jet and a coflowing stream and the resulting

large-scale structures are currently being studied at NBS using a timedependent axisymmetric computer code. When these structures pair and merge,

there is improved mixing between two streams, which is useful in practical

applications. In order to better understand the flow development, controlled

forcing can be applied and the effect of certain types of perturbations on the growth of the mixing layer observed.

Linear stability theory appears useful in predicting some behavior of

the large-scale structures in a jet [1]. The evolution of vortices seems

to be related to the unstable frequency range determined for the upstream

velocity profile. Controlled perturbations applied at the most unstable

frequency and its subharmonics can affect vortex merging and the growth of

the mixing region.

In the present work, the jet's mixing layer is studied using stability

analysis. The objective is to determine the spatial stability of the jet

profile and use the results as the controlled perturbations at the inlet of

the jet computation. Although the procedure developed applies for arbitrary

profiles, the velocity ratio of the jet and freestream has been the main

parameter varied. The analysis assumes that the flow is incompressible,

1
inviscid, and axi symmetric, and that linearized equations are applicable.

The validity of the results can be determined by comparisons of computations and experiments.

In the following section, the equations used in the stability

analysis will be presented. The solution technique using software

currently available on the NBS UNIVAC 1100/82 will then be described.

Results for Gaussian profiles are also included.

Basic Equations

The method presented here parallels the solution procedure described

in reference [2]. The axisymmetric Euler and continuity equations for an incompressible fluid, written in terms of the three variables, u, v, p, are the basis of the formulation, u and v are velocity components in the axial

(z) and radial (r) directions, respectively, in cylindrical coordinates,

and p is the ratio of pressure to constant density. Small disturbances

are added to a known solution: u = U(r) + u' , v = v' , p= p

1 , where the

primed quantities represent the disturbances. Substituting these into

the equations and neglecting higher order terms involving u' and v' yield

the linearized disturbance equations. Assuming parallel flow, the form

of the disturbance is that of a traveling wave whose amplitude varies with

r and which moves parallel to the z-axis.

- *(r)e i(aZ - Bt)

where <j> represents u, v, or p, a is the wavenumber, and 3 , the frequency,

of the disturbances. Substituting this form into the linearized disturbance equations leads to three ordinary differential equations with

complex coefficients.

M - e)v = i ^

(all - B)u - i V ap

For a spatially growing disturbance, a is complex (a + ia- ) and 3 real.

Eliminating u and v from equations (2)- (4) gives an equation for the pressure disturbance.

d£

dr + [- -

L r

2

dU

(u - d

e/a) ]

- a2p = 0

The boundary conditions to be satisfied are that p remains finite at

r = 0 and p becomes zero as r

-* ». Equation (5) and the boundary conditions constitute an eigenvalue problem. The objective is to find the

corresponding values of a and s and associated eigenfunction p(r) for a given profile U(r) and dll(r)/dr.

2
Computational Procedure

Several computer subprograms are used in finding a solution. The

velocity profile being used in the present study is

c

- (r - 1)»

U = 1 - —— e 2 g

z

, r < 1 /2lT

(r - 1)2

u = LI - (IL + — - 1 ) e Zo z 0 0 /27r > 1

where all lengths are nondimensional ized by the jet radius and all velocities by the jet centerline velocity. Equation (6) represents either a single or double Gaussian, depending on the choice of U and C//2tt". U(r)

and dll(r)/dr are continuous at r = 1. The profile U(r)°is specified in the

program as a function UINLET(R) and dll/dr as DUDR(R). Thus, the parameters

or form of equation (6) can be easily changed.

Equation (5) is integrated using the subroutine CDRIV1 from CMLIB

(NBS Core Math Library), which solves initial value problems for complexvalued ordinary differential equations, integrated with respect to a single, real, independent variable. Equation (5) is expressed as a pair

of first-order equations in subroutine F.

yi = y 2

y 2 = [ f

+ 2 ^7 / (U - 3/a)l y 2 + a

2y i

( 6 )

(7)

where y l = p and y 2 = dp/dr. Subroutines USERS, JACOBN, FA, and G, which

are options for the CDRIY package, are included as dummy routines to avoid

error messages during the MAP processing.

Boundary values are needed at r = 0 and at the largest value of r. At each of these locations, dU/dr vanishes and equation (5) becomes a modified Bessel equation. The solutions consistent with the boundary

conditions are p = a I (ar) at r = 0 and p = b K ( a r) as r becomes

large. This leads to Boundary values at r = 0 or p = a 1(0) = a,

dp/dr = act I -. ( 0 ) = 0. Although this result can be determined directly

from axi symmetry, the Bessel function expression is used to obtain the

limiting value of 1/r dp/dr at r = 0 needed for equation (7). Since a is arbitrary, set y^ = 1, y 2 = 0, and integrate from r = 0 to r = 1 to

get the inner solution f^(r). The integration is done in radial increments chosen to provide velocity values at the grid points used in the

jet computation, with the internal step of the integrator adjusted to

just reach each point.

3
The initial values at the largest value of r, r , are p

= bK (or ),

dp/dr = - baK^ar ). Since r is large and evaluating the functions 3 * and with complexarguments is not straightforward, the numerical values

0 are determined from asymptotic expressions. With b also being arbitrary,

set y l = K (ar ), y2 = - a Marmax

)> and inte 9 rate from r = r__„ to

r = 1 , as aescrToed above, to get tne outer solution f (r). max

The required matching conditions at r continuous. This reduces to

= 1 are that p and dp/dr are

df.(D

yu 1dr

d f

0 n)

-Vu-ar- = 0 . ( 8 )

The expression on the left side of equation (8) is written as a complex

function FF whose root, a, is to be found. The integration procedure is

included in the function subprogram.

The values of a are determined by using the subroutine ZANLYT from

IMSL (International Mathematical and Statistical Libraries), which finds zeros of a single nonlinear complex function. For each value of 6 chosen,

ZANLYT is called to find one root, with one guess and no known roots

input. The first guess for a found most appropriate is 6(1 - i/2).

When the root is returned, several tests are made before it is accepted.

The error parameter, IER, must be zero, indicating convergence was obtained. For the problem considered here, other requirements are that

a is positive, a. is negative, and the phase velocity, c , = 6/a , converges to one it 6 = 0. The last condition is checkedp Gy finding a

for a slightly smaller value of 6 (a 6 = 0.01), using the first a as

a guess, and comparing c , values. If c , is less than one and decreases,

or c . is greater than one but less than p T.5 and increases, the value of

6 is p rurther decreased and a determined, to see if the trend reverses in

a short interval (20 A6, A6 = 0.04). If the first root is rejected, a second guess is tried. When convergence was obtained on the first guess

with a positive but a. not negative, the complex conjugate of a is

guessed for the second 1 trial . Otherwise the guess used is 6 + 1/4

+ i / 2 ( 1 - 6) for 6 > 1 or 6(1 - i/10) for 6 < 1. For the cases tested,

a root satisfying the conditions described above was found withoutadditional trials needed.

When determining a as a function of 6, the testing is done only for

the initial value of 6, with A6 = 0.04 in the c , comparison. The guess

for successive points is the value of a from thipreceding 6 value. This

case is of interest for finding the most unstable frequency, the value of

6 for which a- has the largest negative value. The curve is started from

6=4 and proceeds to smaller values of 6, with A6 = 0.04. If the value

of a- is decreasing as 6 decreases, the curve is continued in the same direction until the minimum value of a- has been passed. If increases

in the first set of intervals, the curve is continued similarly in positive increments from 6=4 until the minimum is located.

4
Once a has been found for a particular value of g, the corresponding

velocity eigenfunctions can be determined. The outer pressure solution

f (r) and df (r)/dr are multiplied by the ratio of the constants b/a = f°(l)/f (1) 8 r df. (l)/dr / df (l)/dr, derived from the matching conditions.

Or*ie arbitrary constant remains, multiplying the entire solution. The

velocity v can be computed at each value of r where the integration was

done from the values of dp/dr and equation (2). The u velocity is then

calculated from equation (3) using v and the p solution. In the jet

computation, v is needed at points halfway between those used in the

integration, so averaging is done there. The symmetry conditions at

r = 0 give u(- Ar) = u(Ar), v(0) = 0, v(- Ar) = - v(Ar).

From equation (1), the real velocity disturbance at z = 0 is u‘ =

u (r) cos 3t + u.(r) sin 3t, where u = u + iu. . At a ny value o f r, u

1 i£ a maximum when Bt = tan" 1 (u./u ), giving u' max

= /u| + u. 2

. The

u and u. arrays are searched to determine theradial locations where

each has the largest magnitude. At these two values of r, the quantity

0.1 U(r)/u* is computed and the smaller one selected for the constant

a. The comprete u and v solutions are then multiplied by this factor.

When the procedure is called from the program doing the jet computation, subroutine PERTRB is used. Inputs to this subroutine are the

array of r values where the u velocities are needed and the number of

points in the array, the matching location (r = 1), and S. When 3 = 0

is input, the most unstable frequency is found and that value of 3 is

returned. Additional outputs are a and arrays of u , u . , v , and v •

.

The perturbations applied to the jet velocities are:

N

u' = z a. [u. (r) cos ( 3^t - a. z) + u. (r) sin ( 3^t - k=l

K rK rK

i K

N

v' = z a. [v. (r) cos (3 t t - a. z) + v. (r) sin ( 3^t - k=l

K K r

K K rK

i K

where N frequencies have been included.

Results

Results will now be presented for the inlet velocity profile of

equation (6) for three choices of the constants: (a) a single Gaussian,

U = 0.65, C//27 = 0.35, (b) a double Gaussian, U = 0.3, C I JlH = 0.8,

and (c) a single Gaussian, U = 0.15, C//2tt = 0.85, where a = 0.1 in

each case. These profiles are shown in figure 1. For each of these,

a was determined over a range of 3 * s . For these runs, the stopping

criteria in ZANLYT were (1) the magnitude of the function FF reaching

or falling below 10" 6

, (2) two successive approximations agreeing in

the first six digits, or (3) the number of iterations exceeding the

maximum of 50. The most iterations actually used in finding any one root was 38. The requested relative accuracy in CDRIV1 was set at

10" 6 . The value of r was 5.2.

max

°k z) ^

r

a

k

Z )] ,

5
In figure 2, -a^ is plotted against B> where the variation of the

peak value and most unstable frequency can be seen. The latter ranges

from 3.2 to 4.32 for these cases. In figure 3, c , is plotted against

B. The trend toward c , = 1 at B = 0 is evident nere. The velocity

ratio listed in figures n 2 and 3 is defined as 1/U .

Concl us ion

A collection of subroutines has been developed to solve the

eigenvalue problem resulting from applying stability analysis to the

axi symmetric continuity and Euler equations. The velocity eigenfunctions are directly available for use as perturbations to the

inlet velocity profile in the jet computation.

References

1. Ho, Chih-Ming and Huang, Lein-Saing, Subharmonics and vortex

merqinq in mixing layers, J. Fluid Mech., Vol. 119, p. 443,

1982.

2. Michalke, Alfons and Hermann, Gunter, On the inviscid instability

of a circular jet with external flow, J. Fluid Mech., Vol. 114,

p. 343, 1982.

6
Fig. 1. Velocity profiles
ol/l

"T"

Kj

rvj T

Oo tinQ

t r s s

no_

SpotloL

growth

rot*

vs.

frsqusncy

for

thraa

vstocLly

profits*
Pha»»

velocity

v«.

frequency

for

three

velocity

profiles

NBS-114A (rev. 2-80

U.S. DEPT. OF COMM.

4.

1 BIBLIOGRAPHIC DATA

SHEET (See instructions)

. PUBLICATION OR

REPORT NO.

2. Performing Organ. Report No. 3.

NBSIR 83-2686

Publication Date

May 1985

TITLE AND SUBTITLE

A Numerical Procedure for an Inviscid Stability Analysis of an Axisymmetric Jet

5. AUTHOR(S)

Elizabeth F. Moore

6. PERFORMING ORGANIZATION (If joint or other than NBS, see in struction s) 7. Contract/Grant No. national bureau of standards

DEPARTMENT OF COMMERCE

WASHINGTON, D.C. 20234

8 . Type of Report & Period Covered

Final

9. SPONSORING ORGANIZATION NAME AND COMPLETE ADDRESS (Street. City. State, ZIP)

10.

SUPPLEMENTARY NOTES

| | Document describes a computer program; SF-185, FIPS Software Summary, is attached.

11.

ABSTRACT (A 200-word or less factual summary of most significant information. If document includes a si gnificant

bi bl iography or literature survey, mention it here)

Inviscid stability analysis has been applied to the mixing layer profile of an axisymmetric jet and a coflowing stream. A collection of computer subprograms

has been developed to solve the resulting eigenvalue problem. The effect of

cnanging the velocity profile and its parameters can be easily assessed.

Results for Gaussian profiles are included.

12.

KEY WORDS (Six to twelve entries; alphabetical order; capitalize only proper names; and separate key words by semicolon s) axisymmetric jet; eigenvalue problem; mixing layer; numerical

eigenvalue problems; stability analysis

13. AVAILABILITY

methods for

14. NO. OF

PRINTED PAGES E

SH

Unl imited For Official Distribution. Do Not Release to NT1S

Order From Superintendent of Documents, U.S. Government Printing Office, Washington, D.C. 20402.

11

15. Price

Order From National Technical Information Service (NTIS), Springfieid, VA. 22161 $7.00

USCOMM-DC 6043-P80