A numerical model of acoustic choking, part II: Shocked solutions

Journal of Sound and Vibration (1986) 104(1), 81-107

A N U M E R I C A L M O D E L OF A C O U S T I C C H O K I N G , PART II: S H O C K E D S O L U T I O N S N. J. ~VALKINGTON AND W. EVERSMAN

Department of Mechanical and Aerospace Engineering, University of Missouri-Rolla, Rolla, Missouri 65401, U.S.A. (Received 6 March 1984, and in revised form 4 January 1985) The one dimensional equations of gas dynamics are used to model subsonic acoustic choking. This model can accommodate non-linear distortion of waves and the eventual formation of shock waves. Several finite differencing schemes are adapted to obtain solutions. The results obtained with the various schemes are compared with the asymptotic results available. The results suggest that no one finite differencing scheme gives solutions significantly better than the others and that most of the difference solutions are close to the asymptotic results. If the acoustic shock wave is sufficiently strong it almost annihilates the acoustic wave; in this situation numerical errors may dominate the results. Such solutions involve very large acoustic attenuations.

1. INTRODUCTION In a previous p a p e r [1] the authors considered the non-linear acoustic equations as a special case of the equations of gas dynamics, and used finite difference techniques to approximate their solution. In [1] only continuous solutions were sought~ The finite difference techniques can, in principle, be used to obtain solutions to acoustic problems which involve shock waves. The formation of shocks in acoustic signals can dissipate a significant amount o f their energy and this may have applications in the design of aircraft engine inlets where all of the conditions for significant non-linear behavior are present. In order to design inlets to take advantage of such p h e n o m e n a it will be necessary to develop codes which can handle the non-linear problem with shocks. In this p a p e r we report the resfilts obtained using finite ditterence methods to solve a one dimensional model problem. The presence of shock waves can seriously degrade the accuracy of finite difference solutions. The philosophy normally adopted is that proposed by von N e u m a n n and Richtmyer [2] who introduced the idea of artificial viscosity. Their idea was to accept a degradation in accuracy near the shock wave, i.e., accept a smeared shock, provided that the rest of the solution remains accurate. Since von N e u m a n n and Richtmyer proposed this idea m a n y finite difference schemes have emerged, the motivation being a need to either improve the accuracy, especially near shocks, or to improve efficiency. Presently there is no ideal finite difference scheme. In fact the accuracy of a scheme is usually problem dependent. In this paper the Fubini-Blackstock [3, 4] solution for waves in a uniform duct and the solutions of Myers and Callegari [5-8] and Myers [9] for a TsienCrocco duct are used as standards against which difference solutions a r e compared. Analytical solutions do not, however, provide a completely objective method for comparing the different solutions. Often qualitative comparison of plots is used (see, for example, the p a p e r by Sod [10]). For the problems under consideration here, we compare plots 81 0022-460X/86/010081 +27 $03.00/0 O 1986 Academic Press Inc. (London) Limited

82

N.J. WALK1NGTON AND W. EVERSMAN

and the Fourier coefficients of the acoustic velocity with those of the analytical solutions. Where appropriate the acoustic power is also compared. In the next section we recall the gas dynamic equations o f motion and the b o u n d a r y conditions developed for the acoustics problem in [1]. In section 3 the finite differencing schemes are introduced and in section 4 the results obtained with these schemes are discussed. 2. EQUATIONS OF MOTION AND DUCT GEOMETRY 2.1. D U C T G E O M E T R Y The duct geometry under consideration is shown in Figure 1. For all of the problems discussed here an acoustic source is located at x = 0 and the duct termination at x = l is

Figure 1. Duct geometry.

anechoic (i.e., reflectionless). The ambient flow velocity being negative causes the sound waves to propagate against the mean flow. This is the situation encountered in jet engine inlets where air flows into the engine with the acoustic source modelling the compressor blades and the anechoic termination modelling the atmosphere into which the noise propagates. The negative mean flow tends to " t r a p " the sound in the duct, exacerbating the non-linear effects. 2.2. EQUATIONS OF MOTION The one dimensional gas dynamic equations of motion m a y be written as OU/Ot+OF/Ox=f, where pA

R F=

lp AJ 1EJ

M 2 / R + ~"

,

I(MIR)(E+)J E ---E-,

f=

o}

~rCA'IA) , 0

(1)

in which p, u and e represent the fluid density, v e l o c i t y a n d total energy per unit mass and A the area of the duct. All of the quantities in equation (1) have been reduced to dimensionless forms by using a reference density and speed of sound, taken as the mean flow conditions at x = 0, and a reference length. It is convenient to introduce an alternative spatial co-ordinate, ~:, and an invertible m a p p i n g of ~: onto x to enable difference schemes to be constructed on a uniform sr grid

83

ACOUSTIC C H O K I N G

without any loss of generality. Substituting ~: for x in the equations of motion gives

~V/~t+OF/Or

(2)

where V=(1/Ex)U,F(V,~)=F(U),f=(1/~x)f, andr162 This is the one dimensional form o f the transformation introduced by Vinokur [11]. Finite difference schemes usually require the equations to be modified if they are to yield reasonable approximations to shock waves. Two modifications have been considered. The first simply adds a fourth order term to the equations, giving

OV/Ot + OF/Or =f - v(8r 4 04V/O~4,

(3)

where 8r is the grid spacing on the ~: axis. An alternative to the above is to modify the constitutive model describing the pressure, this being the von Neumann-Richtmyer idea. The modified pressure used here is described by ~" = (3' - 1)(E - M2/R) - ~,,R(lul + a2)(8r 2 0u/0r

(4)

where a is the speed of sound. The constants v and v, depend upon the specific solution being sought and are adjusted to give a satisfactory resolution of the shock. The equations of motion may also be expressed in the non-divergence form,

0 U/at + [A] 0 U/Ox =.7,

(5)

where [A] = [0/3/0 U]. The flux vector, if', is a homogeneous function of degree one; thus it may be expressed as F=[A]U.

(6)

In the lambda scheme introduced in section 4 the following form o f the equations is used:

a lT_I/Ot+ [,,{] a U/Ox = f..

(7)

Here

0

,,

0 =

,

[~] =

0

P = In (p),

0

1

fi=

a:13" u O , it

[ P ] , o

s = In (3"._~)p~,

23"p

=--.p

2.3. BOUNDARY CONDITIONS

The boundary conditions at the duct ends must model the presence of an acoustic source at x = 0 and an anechoic termination at x = l. The boundary conditions may be thought of as specifying that no waves enter the domain 0 < x < l at x = ! and that the waves entering the domain at x = 0 be chosen to model the acoustic source. This is to be accomplished by using the characteristics of the gas dynamic equations. With the left eigenvectors of the coefficient matrix [A] defined in equation (5) defined by /(Or[A] = A(')!c')r,

i = 1, 2, 3,

(8)

the characteristic curves are then defined as those curves in the (x, t) plane satisfying the differential equation

dx(i)/dt = A~i).

(9)

84

N.J.

WALKINGTON

AND

W. E V E R S M A N

The boundary conditions may be summarized as TOU]

= b,

(10)

where [L] = [l m, ! c2~,lt~], and

~l

COT(f-- At~ 0 U/Ox) if the ith characteristic exits from 0 < x < 11 b~ - ~ d a i / d t if the ith characteristic enters 0 < x < l .I The functions ai(t) are set to zero for an anechoic termination and if a source is present a~(t) is chosen to model its presence. This treatment o f the boundaries parallels that used by Hedstrom [12], and has been described in greater detail by Walkington and Eversman [1]. 3. FINITE DIFFERENCE SCHEMES The choice of a finite difference scheme for a particular problem is somewhat arbitrary. There are many schemes to choose from, most o f which have one feature that makes them "better" than the others. One of the first decisions required is the choice between an implicit or explicit scheme. A great advantage of implicit schemes is their ability to take large time steps, i.e., exceed the Courant-Friedrichs-Lewy (CFL) limit. This is very useful for problems which require very fine gnos to resolve sharp spatial features (e.g., a nearly stationary shock or a stagnation point). If such a solution changed slowly the CFL limit would require a time step much smaller than that needed to accurately model its evolution and an explicit scheme may be unnecessarily expensive. However, the acoustics solutions vary rapidly in time as well as space, so a time step similar in magnitude to the spatial grid spacing is required to give uniformly accurate approximations. In this situation explicit schemes are probably more economical. Another decision that must be made is the selection of the order of the scheme. Preliminary experiments with first order schemes indicated that their inherent damping rendered them unsuitable for this problem. Currently second order methods appear to be the most widely used for the gas dynamics equations. Schemes that are of mixed order, usually having greater spatial than temporal accuracy, have been proposed. The advantages of such schemes are similar to those for implicit schemes. If a scheme has a high degree of spatial accuracy it is possible to have a relatively coarse mesh and this will enable larger time steps to be taken without exceeding the CFL limit. However, arguments given below on the use o f hybrid schemes which reduce in order near shocks where derivatives do not exist can be used to rationalize the choice of a second order scheme for a problem where the discontinuous solution is the principal feature. In the present study we have chosen to use proven, robust second order schemes or second order hybrid schemes rather than to experiment with less well tested higher order schemes. Second order schemes are often chosen as candidates for hybridizing in order to try to incorporate within them some logic that enables discontinuities to be modeled better. The schemes of Harten and Zwas [13] and Boris et al. [14-16] are second order schemes with such modifications. Even if a scheme models the smooth portions of the flow field very well the errors introduced by modelling a discontinuity can dominate the whole flow field. This reasoning is used to justify the use of the hybrid schemes which reduce in order near discontinuities. Further justification for such hybridizations can be m a d e by arguing that near discontinuities the solution derivatives do not exist so there is no reason why a high order scheme should be used, especially if it gives poorer answers than a first order scheme.

ACOUSTIC

85

CHOKING

For the acoustics problem a finite difference scheme must satisfy some additional constraints. The scheme must have a steady state (time invariant) solution to the mean flow problem, since the acoustic quantities are defined as the difference between those of the perturbed fluid state and the steady state values. This may appear to be a very mild constraint; however, the MacCormack scheme, for example, is often alternated with the predictor step being forward differenced at one time step and then being backward differenced for the next step. The two versions of the method have slightly different steady state values, which both differ slightly from the exact solution. Alternating them causes the solution to change which results in acoustic noise. Similar problems could be encountered with hybrid schemes which change with the solution. For the acoustics problems under consideration five differencing schemes are considered, all of them explicit second order methods. The schemes selected include the classical differencing scheme introduced by MacCormack [17], the lambda scheme [18, 19] with several variations and the hybrid scheme of Boris and Book (SHASTA) [14]. A second order Godunov method, similar to that proposed by van Leer [20], in which a Lagrangian co-ordinate system is used, was considered by the authors. In this scheme an exact shock tube solution is used which must be solved for each grid point at each time step. This calculation and the necessary remap operations from the Lagrangian to the Eulerian grid resulted in a code that had long run times. The solutions did not appear to be "better" than those of the other schemes so this algorithm was not pursued. Details of this scheme and some of its solutions have been documented in reference [21]. While this is a very small subset of the possible schemes it is expected that it covers a sufficiently large spectrum of those available to show if any differencing method is clearly superior for the acoustics problem under investigation. 3.1. S C H E M E I - - T H E M A C C O R M A C K S C H E M E The MacCormack scheme has been used to solve a wide variety of gas dynamic problems and was used by the authors [1] to solve the continuous problem. In this scheme simple differencing is used to estimate the spatial derivatives:

V; +' = V ; - (StlS~)(F~-+, - F ; ) + 8tf;, v9.~+ , _ ,- ~ ( V j n + V~+,) -(St/28~)(F~+'-F~_*,~)+(St/2)f7

§

(11)

3.2. S C H E M E I I - - T H E L A M B D A S C H E M E Recently Moretti [18] and Moretti and Zannetti [19] introduced the lambda scheme. While this scheme is not conservative they claim that it accurately models steep solutions without oscillations. In this scheme, essentially, the directional derivatives are estimated along each characteristic by using purely one sided differencing and then a linear combination is taken to yield the temporal derivatives of the dependent variables. One defines g~=+(1/48~)(h)

'J•

')+~('),,j•

t'r~ vj -/.'/j• ~) if A~')~ 0,

(12)

where the right-hand side is evaluated at the nth time step, and then evaluates the half step values as [L];r[/]~+~- 0~ + (St/2)f~.] = ( - S t / 2 ) ~ x , g ~ .

(13)

Defining g~§ similarly by using the half step values and then putting

G,j-g,j2+(g,j-gij• __

n+

!

n

n

ifA)')~ 0

(14)

86

N . J . ~,VALKINGTON A N D W. E V E R S M A N

enables the full step to be written as [ [.]~,+,T( 0~,+1 _ 0~' +

8tf; +~)= -St#=,Gj.

(]5)

3 . 3 . S C H E M E I I 1 - - T H E SPLIT FLUX S C H E M E

Recently Steger and Warming [22] introduced a conservative version of the lambda scheme. Their idea was to split the flux vector, F, into two components, F + and F - , each of which can be associated with the characteristics having positive and negative slopes respectively. This is done by decomposing the coefficient matrix [A] in equation (5) into the sum of two matrices each of which has either positive or negative eigenvalues. Explicit formulae for F + and F - have been given in reference [22]. The difference scheme is then written as V ; + t = V ; - ( St l S f )[ F;+ , - F ; + F ; - F+_,] + ~tJ ; ,

V~ '+' = ~(V~' + V ; + ' ) - ( & t / 2 ) f ;

-(F;§

+' - ( S t / 4 $ ~ ) { ( F ; + , - F ; ) o + I + ( F f + , - F f ) "

-F;.+,)"+'+(F;--F;+-,)"-(F~+-,-F;_2)"}.

(]6)

3.4. SCHEME IV--A MIXED SCHEME

Steger and Warming [22] found that alternating schemes I and III with each time step produced an algorithm that appeared to be better than each individual scheme. This is not possible for the acoustics problem since a steady mean flow is not attained; however an analogous hybrid scheme may be constructed simply by averaging the two schemes: i.e., vylltv

-_ ~,,(, /v.,,+,l , , + V~'+'lm),

(17)

where the sub index indicates the scheme used. 3. 5. S C H E M E V

T H E SHASTA A L G O R I T H M

Boris, Book et al. [14-16] introduced the SHASTA scheme. This scheme was designed to eliminate the spurious oscillations that often occur near shock waves. To accomplish this the scheme reduces to order one whenever certain monotonicity constraints are violated. The SHASTA algorithm is somewhat involved; details can be found in reference [14]. 3.6. v i s c o u s TERMS The fourth order viscous term described in equation (3) is incorporated into the above differencing schemes according to V~+llncw=

V~+l[ol~l-v[ V;+2-4 V~+l+6 V~-4 V~-i + V]-2].

(18)

When implementing the von Neumann Richtmyer type viscosity the derivative in equation (4) was estimated by using a central difference. The viscous terms were not applied at points on or adjacent to the boundaries. This avoids the practical problem of having to estimate the values of the dependent variables at points which would lie outside of the domain 0 ~< x<~ 1 as well as the theoretical problem associated with the application of the correct boundary conditions when the equations are no longer hyperbolic. 3.7. BOUNDARY POINTS The differencing schemes I to V may be used for points in the interior of the domain; at the boundaries different algorithms which incorporate the boundary conditions are required. The values of the dependent variables on the boundaries were found by using

87

ACOUSTIC CHOKING

a Taylor's series, u~ § = u~ + 0 u/atl~ ~t + ~2u/ot2l~

8t~/2.

(19)

The first temporal derivative may be evaluated directly from equation (10) and the second temporal derivative obtained by differentiating equation (10) with respect to time. Schemes II and III need to be modified for points adjacent to the boundaries since they require two neighboring points. For the first interior point scheme III requires the + 4term (Fo - F - 0 to be estimated. When expanding about point one the term becomes ( F ~ - F+_O = aF+la~l • 8~-~(o2F+1c3~ 2) 8~ 2 = 2aF4-/a~[~ 8 ~ - (F~ - F~).

(20)

The derivative in the latter expression may be evaluated from the equation of motion, aF§

= f - a V/at - aF-/as

(21)

A similar procedure is used for the boundary at x = l and to evaluate the g~. terms for the lambda scheme if the corresponding characteristic enters the domain. 3.8.

SOLUTION

PROCEDURE

The algorithms discussed in the previous section will evaluate the transient response to a signal travelling into a duct; however almost all of the analytical results are for periodic solutions. Accordingly the following procedure was used to obtain similar results. (a) Use a Runge Kutta algorithm to generate an approximate mean flow profile for the initial conditions. (b) Run the program with zero boundary conditions until the solution changes negligiblyin order to obtain the steady state solution. (c) Introduce the periodic acoustic signal, via the boundary conditions, and sustain it until all the transients propagate out of the ends of the duct. Store the periodic solution. (d) By using a fast Fourier transform (FFT) routine the temporal solution at any position is decomposed into its harmonic components. Two temporal samples, one quarter o f a period apart were used to check that the solution was periodic. One problem with this approach is the possibility of a drift in the steady mean flow established in step (b). It is possible that an analytical mean flow be used and the non-linear acoustic equations be solved numerically with the analytical mean flow appearing as specified coefficients. We chose not to follow this approach for two reasons. First, the analytical solutions are available only for one dimensional flows. The one dimensional model discussed here is intended to investigate the viability o f solving non-linear acoustic problems, the ultimate goal being the solution of the multi-dimensional problem. Problems associated with the specification o f a mean flow are likely to arise in the multi-dimensional problem and the one dimensional problem yields information as to what they are. Second, a shift in the mean flow is a solution to the non-linear acoustic equations and their finite difference counterparts. Since both the periodic solution and the steady state solution are approached from a transient finite difference solution it is important that the steady state solution be absolutely consistent with the periodic solution with respect to the mean flow component or errors of acoustic magnitude may appear. 4. RESULTS The Fourier components of the acoustic velocity at x = I are reported here. The acoustic velocity at the anechoic duct termination may be reconstructed according to the equation

88

N.J.

WALKINGTON

AND

W. E V E R S M A N

(This is the signal which would exit from a jet engine inlet into the environment.) When appropriate, the acoustic power at the duct terminations will be reported. For all of the results presented the ratio o f specific heats is taken as y = 1-4. When harmonic amplitudes and frequencies are converted to sound pressure levels (SPL's) and dimensional frequencies ( f ) a reference density of pref = 1"2 kg m -3, a speed o f sound of Crer= 340 m s -I and length Lrer= 1 m are used. The sound pressure levels and frequencies then correspond approximately to those in air at sea level. 4.1.

RESULTS FOR A UNIFORM

DUCT

A uniform duct o f length ! = 0.5 with a source amplitude a I = 0.02, frequency o = 2~r, and an adverse mean flow Mach number of M = - 0 . 7 5 was selected for the first example, Case 1. This combination of parameters corresponds approximately to a sound pressure level SPL = 160 dB and frequency f = 340 Hz and results in a shock forming at x = 0.41. The shock wave which forms propagates in the positive x direction and exits the duct at x = 0.5. This problem will then require the boundary conditions, originally developed for continuous solutions, to pass a shock wave. The solution to this problem becomes (infinitely) steep towards the right-hand end of the duct, so a mesh transformation was used to concentrate the grid points toward this end. A quadratic transformation of the form

x/I=[(2X+I)-x~?]/(X+I),

r/= (~:+ 1)/2,

-1<~:<~ 1,

(23)

was used to accomplish this. The parameter X alters the degree o f distortion. If )t' = 0 the mesh is uniform and for X > 0 the grid points concentrate at the right hand end. For this problem X was set to ~. The resulting grid has grid spacing at x = ! of about 0.43 of that at x = 0 . The time step was chosen so that 8t/SXlmax=l and a mesh of 201 grid points was used. The harmonic velocity amplitudes at x = 0.5 for each of the schemes I-V with zero viscous terms are shown in Table 1. Included in Table 1 are the Fourier coefficients TABLE 1

Fourier coefficients for case 1 with zero viscous terms ( N N =201, X = 2 / 3 ) xlO 2

Blackstock'~ Scheme I

.

Ib.I

Ib.I

0~

1 23 4 5 6

1"659 0"721 0"449 0.322 0.250 0.204

1-650 0-712 0.440 0-316 0-247 0-207

-92 -97 -100 -102 -104 -106

Scheme II

Scheme III

Scheme IV

Scheme V

Ib.I

Ib.I

Ib.I

Ib.I

1"668 0.734 0"466 0.346 0"280 0.238

o: -91 -93 -93 -91 -87 -90

1.672 0.740 0-471 0.345 0.269 0-215

0:~ -91 -94 -93 -91 -87 -87

1"663 0.730 0"461 0-337 0-267 0.220

07, -92 -96 -97 -98 -98 -99

1.651 0.713 0.436 0"306 0.234 0"189

0~, -92 -96 -98 -99 -101 -104

t 0. = - 9 0 ~ f o r all n.

calgulated by using the Fubini-Blackstock solution. Plots o f the acoustic signal emanating from the duct at x = 0.5 over one period are shown in Figure 2 for each scheme. Plotted with the difference results are the linear (sinusoidal) and Fubini-Blackstock solutions. (The Fubini-Blackstock solution was calculated by using the parametric equation that comes from the characteristic solution since the corresponding Fourier series converges very slowly near the shock and the calculation of the coefficients is not trivial.) It is apparent from Figure 2 that scheme I and to a lesser extent scheme II do not yield good

ACOUSTIC

89

CHOKING

o >.

~176 / \

9 0 020

O010

0 030

t /

"\ "~',

.o o 0 1 0

4 25 4 50k~'~.. 4"~'5 5, tOO

0 000 -0010

(o) \

\

~

-0 020 0020 0015 :~" 0010 0005 _'2 o oo0

0 020 L

~-'~ " ~

V

0 020 ,

"\

0015 0 010 0005

o ooo 0 005

~-0005

~-0010 -:0015

4 25

0010

-0 020 -0025 -0 030

0015

(b) "N.~'

0-020

(d)

~"

OO2O O0 15 -->" 0010

9

o

._u

005

0 000

8 -0005

450%,

T' e'xXX

5 00

/

-0010 -0 015 -0 020

Figure 2. Acoustic velocity at x = I for case 1. N N = 201, X = 2 / 3 . ( a ) S c h e m e I; ( b ) scheme II; (c) scheme 1II; (d) scheme IV; (e) scheme V. - - - , F u b i n i - B l a c k s t o c k ; - - - - - , linear; , numerical solution.

approximations to the shock wave. This is a situation where one must "adjust" the viscous terms if the schemes are to "satisfactorily" approximate the solution. Figure 3 shows the solutions obtained by using non-zero viscous coefficients with the differencing schemes. Figure 3(a) shows the solution obtained by using scheme I and a von Neumann-Richtmyer type viscosity with a viscous coefficient of v, = 0.2. Figures 3(b)-(e) show the solutions obtained by using schemes I-IV with the fourth order viscous term and l,= 0.02. The viscous terms were not used with scheme V since it was developed primarily to eliminate the need to use such terms. The harmonic amplitudes obtained for these solutions are shown in Table 2. It is clear from the harmonic amplitudes and the plots that this particular choice of fourth order term changes the solutions obtained with schemes III and IV negligibly. However, the fourth order term made a striking improvement in the solution given by scheme I and a mild improvement in the solution given by scheme II. It is apparent from Figure 3(a) that the von Neumann-Richtmyer viscous term tended to reduce the amplitude of the spurious oscillations but did not reduce their extent like the fourth order term does. Accordingly only the fourth order term is used for the subsequent examples.

90

N. J. WALKINGTON AND W. EVERSMAN

oo3oliI

o )2C

0 025~-[~

0 )I[

oo2obl-l-#'k..,~.

_>.oo,II

~ OOlOI-t/

o| oooo~f

0 )IC 0 )0s

\,

0 )OC ,

g | ta -O 0051< 20010 [-

4-25

-oo,51-

"%

,

4 50 ,w..,,~ Timo \ ~ .... \,,'~

4.75 /

x.~i

,_, \.

_0 ozoL

"... j -

lul

'

-,-4.25

4, 5 50o

~OC '0 )0[ , "0 )IC ,

/

"~..gs

M

\\

",,\,i {

0 )I (c)

0 )2C

, ~ I ~ '

OOEC

0025 00EO ~, 0 0 1 5

o01~

0010 I ~> 0 0 0 5

0005

001C

//

~

,

,

:~ 0.000

O00C

4~5

,

,50kx"

475

5001

0005

~-ooo5 I -0010

0010

-0015

0015

0 0201

0 020

(b)

o o15~

,/...,~

.,=_ 0.010'

oooo

,. ~,

4.~5 450",~,,

4%

5,;r

-0010 -0015 (e)

-0020

Figure 3. Acoustic velocity at x = I for case 1 with non-zero viscous terms. N N = 201, X = 2/3. (a) Scheme I, v . = 0 . 2 0 ; (b) scheme I, t,=0.02; (d) scheme II, v=0.02; (d) scheme III, v=0.02; - - - , Fubini; linear; , numerical solution.

TABLE 2

Fourier coefficients f o r case 1 with non-zero viscous terms ( N N = 201, X = 2/3) Scheme I v , = 0.20~

Scheme I v = 0.02

S c h e m e II v = 0.02

S c h e m e III v = 0"02

S c h e m e IV v = 0.02

[b.I

]b.I

o.

lb.[

o~

Ib.I

o:

lb.I

o~

Ib.I

o~

1"659 0"721 0"449 0"322 0.250 0"204

1"640 0.706 0.434 0.310 0.241 0"199

-92 -97 -100 -103 -105 -109

-92 -97 -100 -103 -106 -111

1.672 0.738 0.469 0.346 0.273 0"222

-91 -93 -92 -88 -83 -77

1"674 0.742 0"470 0"340 0"259 0"200

-91 -93 -92 -89 -86 -81

1.666 0-733 0.463 0.338 0.264 0.212

-92 -95 -96 -97 -97 -97

xIO 2

Blackstockt

n

1 2 3 4 5 6

t O. = - 9 0 ~ for all n.

1-656 0"723 0.454 0.331 0.262 0"218

ACOUSTIC

91

CHOKING

It is apparent that all the finite differencing schemes, with the appropriate viscous terms, given reasonable approximations to th~ acoustics problem and that no one shines over the other. The boundary condition passed the shock wave without any apparent problems. This fortuitous result may be due to the fact that the numerical algorithms all represent the discontinuity with a steep but continuous solution, which may enable boundary conditions developed for continuous solutions to suffice. In all of the problems considered there was never any sign of spurious reflections at the duct terminations even when a shock wave was at the exit. Figure 4(a) shows the spatial variation o f acoustic velocity given by scheme I ( v = 0 . 0 2 ) just as the shock is at the right-hand end of the duct (t = 5.0). Figure 4(b) shows the spatial variation o f acoustic velocity given by scheme V at the same time. This figure shows that while scheme V may eliminate spurious extrema from forming near the shock wave it does not discriminate between spurious and valid extrema, so all extrema are incorrectly modeled.

0015

.oo,oo, _; ,/,'t , H 0~00/ m0 0 0 5

0200~ 0300 J0400

/

Distance I~oIong

5oo

/]dud

-~176176 U -0015

-0 020

(o }

o ozo

/,,~%

//

o <-ooo5 -ooioi-\

, /

\

\/

/ 0000

//\

'

O.OLO

9

,,-,/-

//\\

=,. 0 0 1 5

'

i/

y I

,

I

,

O.lOO1 o-zoo~ 0300 io4oo /, o,s,o~e l\o,o~g//d~c, It

I',

o'5oo ,

I:

F i g u r e 4. S p a t i a l v a r i a t i o n o f a c o u s t i c v e l o c i t y f o r c a s e 1 at t = 5.0. N N = 2 0 1 , X = 2 / 3 . (a) S c h e m e I, v = 0 . 0 2 ; ( b ) s c h e m e V. - - - , F u b i n i - B l a c k s t o c k ; - - - - - , linear; ---, numerical solution.

4.2. N O N - U N I F O R M D U C T S O L U T I O N S The non-uniform ducts are chosen to be converging-diverging nozzles. Near the throat the mean flow Mach number can be made very close to unity; the non-linear effects will then be significant in this region. At the two ends of the duct the mean flow is usually significantly less than unity and the harmonic amplitudes are small. In this situation it is reasonable to use linear acoustic theory for acoustic power calculations at the duct terminations. Because the solution changes very rapidly near the throat a mesh transformation was chosen that would concentrate the nodes in this region. A cubic transformation of the form

x/l=89162162

- 1 ~< ~:~<1,

(24)

92

N . J . W A L K I N G T O N A N D W. EVERSMAN

will give a suitable mesh. The free parameter X plays a similar role to the parameter X in equation (23). When X = 0 the mesh is uniform and for all positive values o f X the grid points are concentrated near the throat. For all o f the examples X was chosen as either X =89 or X = 1.0. The ratio of grid spacing at the throat to that at the duct exit then becomes

8x,18 oI = , --89 or For a fixed n u m b e r o f grid points and fixed C F L n u m b e r the number of time steps requ!red for X = 1.0 is approximately one and a half that required for X = 1/3 since the finer grid near the throat necessitates a smaller step size. The non-uniform duct solutions are considerably more expensive (computationally) than the uniform duct solution. For all of the non-uniform results presented the steady flow was given 20 (dimensionless) seconds to establish. The Mach number at the throat was printed out every two seconds and after 20 seconds it was observed that the first five (and usually six) decimal digits did not change. The periodic signal was then introduced at the left-hand b o u n d a r y and sustained for 30 seconds. The temporal solution in the final periods was stored and then post processed with the FFT routine. Two temporal samples approximately one quarter o f a period apart were processed with the F F T routine. This procedure produced Fourier coefficients for each sample whose fourth significant figure differed at most by two. Preliminary experiments indicated that at times much less than 30 seconds the solution was not periodic since the Fourier coefficients for the two samples differed significantly. 4.3. QUARTIC D U C T I f the source amplitude and frequency are small and the Mach number is significantly different from unity the non-linear solution should almost exactly coincide with the linear solution. The non-uniform duct solutions involve a non-zero non-homogeneous term and the acoustic signal interacts with the flow in a n o n - t r i v i a l way. In order to verify that each scheme would model these p h e n o m e n a an almost linear case was run for a duct with a quartic area variation given by

A(x)= 89

0<~x~
(25)

The inlet/outlet Mach n u m b e r was chosen to be M = - 0 . 3 which results in a throat Mach number o f M = -0.86. A velocity amplitude of a~ = 2.0E - 0 5 (100 dB) at a frequency of to = 1.9 ( f = 103 Hz) were selected for the source at x--- 0, and a mesh transformation parameter o f x =~ on a mesh with 151 grid points was chosen. The m a x i m u m ratio 6t/6x was approximately 0.4. Table 3 shows the transmitted harmonic amplitude for the first harmonic for each scheme. The higher harmonics were all very small. Also included in Table 3 are the C P U times (Amdahl 470; included here only for comparison of schemes) and the steady throat Mach number given by each scheme. It is apparent that schemes I - I V all g i v e very similar answers while scheme V gives slightly larger deviations from the linear solution and exact throat Mach number. The reason for this discrepancy becomes apparent when one views the plot of the spatial variation o f acoustic velocity for a fixed time. Figure 5 shows such a plot at t = 30 for scheme V with the linear solution plotted on top. It is clear that the extrema are modeled very poorly. It was observed that this problem originates from the interaction of the extrema in the acoustics solution with the extrema in the mean flow. For problems that have highe r frequencies, and hence more extrema that are close together, scheme V would not give a plausible solution, van Albada et al. [23] reported a similar failure of scheme V. While schemes I - I V all gave similar

93

ACOUSTIC CHOKING TABLE 3

Results for case 2 ( N N = 151, X = 1/3) Transmitted amplitude Scheme

CPU (min)

Throat Mach number

Ib,lx 10'

o;

I II III IV V Linear

5-4 8-3 9.2 11.8 14.7 --

-0-862 -0.859 -0.859 -0.860 -0-855 -0-861t

1-540 1-539 1-537 1.538 1-511 1-539

-2.0 -0.5 -0.5 -1.2 2.8 -1.8

^

t Exact throat Mach number.

xlO

"4

1.500 1-250 1.000 0 -750 0 -500 0-250 .9 0-000 -0-250 -0-500 b

0

o25~

I

o50/

,r

Distonce~~duct

I

o-r5

I

1.oc

-O'750 - 1000

-1.250 - 1-500, Figure 5. Acoustic velocity profile given by scheme V for case 2.

results the C P U times are not at all similar. Schemes I and II were considerably cheaper than schemes IV and V. 4.4. TSIEN-CROCCO DUCT Myers and Callegari [5-8] and Myers [9] developed non-linear solutions for a Tsien Crocco duct using the method of matched asymptotic expansions. The Tsien-Crocco duct profile has a linear m e a n flow Mach number when the throat is just choked. An analytical solution for the linear acoustic field is available for this configuration expressible as a linear combination of hypergeometric functions. A description o f the mean flow and the acoustic solution may be found in Appendix 2 o f reference [6]. All of the results reported are for a duct of unit length and have an inlet to throat area ratio of 1.32. Figure 6 shows the duct area as a function of position. Four different solutions are presented for this duct profile. The difference solutions are c o m p a r e d with asymptotic results to verify their integrity. : All o f the results for the Tsien-Crocco duct have the velocity source chosen so that the periodic acoustic velocity at x = 0 is a pure sinusoid. This requires the source to have higher frequency components to "cancel" the corresponding reflected waves. However, the non-uniform duct profiles with a high mean flow Mach n u m b e r produce almost no reflections so the superharmonic components required in the source are very small. The

94

N. J. W A L K I N G T O N A N D x,V. E V E R S M A N 15 0

J

1-25 100 0.75 0-50 0 25

0 00

025

050 )

lOO

0-75

Distance along duct

Figure 6. Area variation for the Tsien-Crocco duct. TABLE 4

Acoustic source and mean flow data for the Tsien-Crocco duct solutions (cases 3 - 6 ) Throat Mach Case number 3 4 5 6

-0.975 -0"95 -0.95 -0.98

, Amplitude 2.435 7"385 2"377 2.421

Velocity sourc~ " . SPL(dB) Frequency to f ( H z )

E-4 E-4 E-3 E-3

122 131 142 142

1.026 2"051 2"051 2"051

55 110 110 110

Propagating power " , x=0 x=l Ratio 9.42 8"65 8"98 9"20

E-9 E-8 E-7 E-7

9.42 E - 9 8.44 E - 8 2-57 E - 8 1"03 E - 8

1-00 0.98 0.29 0"01

m e a n flow a n d s o u r c e d a t a for cases ( 3 - 6 ) a r e s h o w n in T a b l e 4. T h e ratio 6t/6Xlmax was m a i n t a i n e d as 0.4 f o r all o f the T s i e n - C r o c c o d u c t s o l u t i o n s . S c h e m e V was n o t u s e d for these p r o b l e m s . T h e s o l u t i o n to c a s e 3 is highly n o n - l i n e a r b u t c o n t i n u o u s . It is i m p o s s i b l e to tell that the w a v e h a s not s h o c k e d s i m p l y b y l o o k i n g at a p l o t o f the solution. S c h e m e I gave s o l u t i o n s to this p r o b l e m that h a d o s c i l l a t i o n s n e a r the steep p o r t i o n o f the wave. A v i s c o u s t e r m o f v = 0.01 was u s e d to e l i m i n a t e this. S c h e m e s I I - I V all have sufficient i n h e r e n t d a m p i n g to give n o n - o s c i l l a t o r y solutions. T h e t r a n s m i t t e d h a r m o n i c a m p l i t u d e s o b t a i n e d with t h e differencing s c h e m e s for case 3 on a m e s h with 251 grid p o i n t s a n d X = 1/3 are t a b u l a t e d in T a b l e 5. T a b l e 6 lists the F o u r i e r coefficients o b t a i n e d with s c h e m e s 1 - I V on a m e s h with 301 nodes. C o m p a r i s o n o f the results in T a b l e 6 with those in T a b l e 5 shows h o w t h e finer mesh is b e t t e r a b l e to TABLE 5

Fourier coefficients for case 3--coa~e mesh (X = 1/3) Scheme I

N N =251 xl04 n

1 2 3 4 5 6 7

Myers

Ib.l s

2.232 0"770 0.405 0.253 0"176 0"134 0"107

v=0.01

Scheme II NN=251

Scheme III NN=251

Scheme IV NN=251

0;

Ibd

0;

Ib.I

0~

Ib.I

0~

Ibd

0;

100 -80 114 -50 149 -11 -173

2-233 0"723 0.344 0.180 0"096 0.050 0"026

100 -79 112 -60 127 -49 134

2.253 0"948 0.372 0.209 0.123 0"075 0"046

101 -74 127 -28 -178 35. -112

2-200 0"581 0.185 0.058 0.017 0"005 0.002

103 -75 116 -54 136 -34 156

2.222 0"655 0"253 0"099 0.037 0.075 0.046

102 -77 144 -58 130 -42 146

95

ACOUSTIC C H O K I N G TABLE 6

Fourier coefficien ts for case 3--fine mesh ( N N = 301)

xlO 4

n 1 2 3 4 5 6 7

Myers

Ib.I

0:

2-232 100 0.770 - 8 0 0.405 114 0-253 - 5 0 0-176 149 0-134 - 1 6 0.107 -173

Scheme I X = 1/3 v=O.O1

0.

x=l]3 lb.I o~

101 -79 114 -54 136 -36 149

2-249 101 0-752 -75 0-382 124 0"222 -32 0"137 176 0"089 26 0.059 -123

Ib.I 2-236 0"736 0.363 0.203 0.120 0.072 0.042

Scheme II

Scheme III x=l/3

Ib.I

0:

2-220 102 0.640 - 7 5 0-238 116 0.089 - 5 3 0"032 137 0"011 -33 0"004 157

Scheme IV

X=I/3 Ib.l 0~, 2-232 0-692 0"298 0.135 0.060 0.026 0"011

102 -77 115 -54 135 -36 153

Scheme 1 X = 1.0 v=O.Ol

Scheme ]I x=l.O

Ib.I

0:

Ib.I

2.239 0.749 0"385 0"230 0"151 0"106 0"076

101 -79 115 -50 148 -18 173

0:

2.244 101 0"755 - 7 5 0"391 120 0-237 - 3 8 0.157 167 0"111 14 0"082 -139

resolve the higher harmonics. Also listed in Table 6 are the results obtained with schemes I and II on a mesh of 301 nodes with a mesh transformation coefficient X = 1.0. This coarsens the mesh at the duct inlet and outlet but refines the mesh in the throat. The subsequent improvement in accuracy of the results suggests that the resolution in the throat is limiting the accuracy of the solution. Only schemes I and II could be run with this mesh since the other schemes would require more than an hour of processor time with this configuration and this was beyond the resources available. Plots of the acoustic velocity obtained by schemes I and II on the fine mesh are given in Figures 7 and 8. The plots of acoustic velocity as a function of position show how large the solution becomes in the throat when it encounters the high Mach number flow. It would be difficult to devise a shock detection algorithm to be used for shock fitting that does not detect a shock in the solutions shown in Figures 7 and 8. x l O -3 9 60001 5-00014000]>. 3-O00J"~

210 0 0 ~

u

0 O00V

(a)

~oool- -I

025

-000}-

-2ooo i-

1.'oo

0 Distance

l 310 0 0 ~ l

-4-000~-5-000~-6000 L

xlO -4 2"000 1-500 1-000

. ,\,

~;~\, (b)

9 o-sool o 0000 --O "5 0 C

8 <{

z3oo

I

,~5oo \k

I

2700, Time

I

2900

/

- I -ooo -1.500 -2-000 -2" 500

,

/

Figure 7. Solution to case 3 given by scheme I. NN = 301, X = 1.0, v = 0.01. (a) Spatial variation; (b) temporal variation at x - / . - - -, Myers; - - - - - , linear; - - , numerical solution.

N. J. WALKINGTON AND W. EVERSMAN

96

xlO -s 6 000 5 000 4 000

/

>, 3 000 z ooo ~ooo

..=

u 0 000 - ~ ooo ~-2000

(a)

i

i

0 25

0 5C Distance

-3 000 -4 000 -5 0 0 0

t

1oo

F

long duct

-6 000

1-500 ]

o 0"~176176 500 -~ 0 000 -o 500 r

obo -1"500 r -2 000~

-2500[Figure 8. Solution to case 3 given by scheme II. NN = 301, g = 1.0, u = 0-0. (a) Spatial variation; (b) temporal variation at x =/. - - -, Myers; - - - --, linear; , numerical solutions.

T h e acoustic p o w e r at x = ! was calculated for each o f the solutions given in Tables 5 a n d 6. The powers for the coarse mesh are listed in T a b l e 7 a n d those for the fine mesh in T a b l e 8. T a b l e s 7 a n d 8 show that all o f the differencing schemes give less total p o w e r t r a n s m i t t e d t h a n the a s y m p t o t i c results o f Myers. H o w e v e r as the mesh is refined the difference results get close to the a s y m p t o t i c results. The powers were calculated a c c o r d i n g to the e q u a t i o n ~w- = ~ 1 ,

I

Re

(26)

tpoao where u, are the h a r m o n i c velocity a m p l i t u d e s (b, at x = 1), p , are the h a r m o n i c pressure a m p l i t u d e s , a n d the subscript 0 indicates a m e a n flow q u a n t i t y .

TABLE 7

Transmitted powers for case 3--coarse mesh (X = 1/3)

n

Myers

NN=251, v=0.01

Scheme II NN=251

Scheme III N N - - 251

Scheme IV NN=251

1 2 3 4 5 6 Total

7"96 0-95 0.26 0.10 0"05 0"03 9.42

7-97 0"83 0.19 0"05 0.01 0.00 9.06

8" 11 0"89 0.22 0-07 0.02 0"01 9-33

7.74 0.54 0.05 0.0l 0.00 0.00 8.34

7.89 0-68 0.10 0.02 0"00 0"00 8-70

Scheme I

xl0 9

97

ACOUSTIC CHOKING TABLE 8

Transmitted powers for case 3--fine mesh ( N N = 301) Scheme I •

X=]

Schemell

n

Myers

~, =0.01

X =t

SchemelIl X =I

1 2 3 4 5 6 Total

7-96 0-95 0.26 0.10 0.05 0.03 9-42

7.99 0.86 0.21 0.07 0.02 0.01 9.16

8-08 0.90 0.23 0.08 0.03 0.01 9-35

7.88 0-66 0.09 0-01 0.00 0.00 8.64

Scheme IV X =~

Scheme I x=l-0 v =0.01

SchemelI X = 1.0

7.96 0.77 0.14 0.03 0.01 0.00 8.90

8.01 0.90 0.24 0.08 0.04 0.02 9.30

8.05 0.91 0-24 0.09 0.04 0.02 9.37

T h e results for cases 1-3 d o not s h o w a n y o n e o f the d i f f e r e n c i n g schemes to be s u p e r i o r to the others. S c h e m e s I a n d II require less p r o c e s s o r time f o r a given grid a n d d o n o t give results i n f e r i o r to the others. The results also i n d i c a t e t h a t it is necessary to use a fine grid in o r d e r to resolve the h i g h e r h a r m o n i c s no m a t t e r w h a t s c h e m e is used. S c h e m e s I a n d II a r e t h e n m o r e attractive since for a fixed p r o c e s s i n g time they can be u s e d o n finer grids. In o r d e r to m o d e r a t e the c o m p u t a t i o n a l e x p e n s e o n l y s c h e m e s I a n d II w e r e u s e d for cases 4 : 6 . F o r all o f the s h o c k e d results scheme I r e q u i r e d a viscous term. S c h e m e II r e q u i r e d a viscous t e r m o n l y for cases 5 a n d 6, these two cases h a v i n g s t r o n g e r shocks. G e n e r a l l y it was f o u n d t h a t the m o s t a c c u r a t e results were o b t a i n e d with the m i n i m u m p o s s i b l e viscous t e r m that w o u l d localize the n o n - p h y s i c a l o s c i l l a t i o n s a r o u n d the s h o c k wave. W h e n t h e viscous t e r m s were not used f o r p r o b l e m s i n v o l v i n g strong shock w a v e s t h e s o l u t i o n s were n o t p e r i o d i c . I f m o r e t h a n the m i n i m u m v i s c o u s coefficient was u s e d t h e s h o c k e d s o l u t i o n s g e n e r a l l y gave h a r m o n i c a m p l i t u d e s that w e r e t o o large for the l o w e r h a r m o n i c s a n d t o o s m a l l for the higher h a r m o n i c s . T h e l a r g e r a m p l i t u d e s o f t h e l o w e r h a r m o n i c s c a u s e d the t r a n s m i t t e d p o w e r to be o v e r e s t i m a t e d . As the mesh was refined t h e t r a n s m i t t e d p o w e r p r e d i c t e d b y the n u m e r i c a l s c h e m e s d e c r e a s e d t o w a r d s that given b y the a s y m p t o t i c results. This is in c o n t r a s t to the results o b t a i n e d for c o n t i n u o u s s o l u t i o n s w h e r e the viscous t e r m s t e n d to reduce t h e p o w e r t r a n s m i t t e d . This suggests that t h e TABLE 9

Fourier coejficients for case 4--coarse mesh ( N N = 251, X = 1/3)

X 10 4

Myers

Scheme I v = 0-01

Scheme II ~,= 0.00

Scheme II v = 0-01

.

Ibol

o:

Ib.I

o•

Ib.I

o:

Ibol

o:

1 2 3 4 5 6 7 8 9 10

6.272 2"579 1.540 1"083 0"832 0.675 0"568 0-490 0-431 0.347

179 80 -13 -103 166 76 -15 -105 164 73

6-351 2.576 1"520 1"033 0.718 0"484 0"313 0"197 0"121 0"074

-179 83 -10 -104 151 56 -48 -152 104 -1

6.405 2.628 1-563 1.058 0.704 0"511 0-342 0.219 0"135 0.081

-178 90 9 -65 -135 159 96 34 -28 -89

6.419 2"596 1.463 0-895 0-546 0"322 0.184 0-101 0.055 0"029

-178 91 11 -64 -137 152 82 13 -56 -123

N. J. WALKINGTON AND ~,V. EVERSMAN

98

TABLE 10

Transmitted powers for case 4--coarse mesh ( N N = 251, X = 1 / 3 ) xlO s n

Myers

Scheme I v = 0-01

S c h e m e II v = 0.00

Scheme II v = 0.01

1 2 3 4 5 6 7 8 9 10 Total

6-32 1"07 0.38 0"19 0.11 0.07 0.05 0"04 0.03 0-02 8.44

6"48 1.07 0"37 0"17 0.08 0"04 0"02 0.01 0"00 0"00 8.23

6-59 1.11 0.39 0.18 0.09 0.04 0.02 0.01 0-00 0-00 8.43

6.62 1"08 0-34 0.13 0"05 0"02 0-01 0"00 0"00 0"00 8.24

TABLE 11

Fourier coefficients for case 4--fine mesh ( N N = 3 0 1 )

xlO a

Myers

Scheme I v=O'O1 X = 1/3

S c h e m e II

v=O'O X = 1/3

S c h e m e II v=O'O1 X = 1/3

Scheme I v=O'Ol X = 1.0

S c h e m e II v=O'O X = 1-0

n

Ib.I

o:,

Ib.I

o:,

Ib.I

o~,

lb.I

o~,

Ib.I

o:

Ib.I

o:

1 2 3 4 5 6 7 8 9 10

6.272 2.579 1-540 1.083 0.832 0.675 0.568 0.490 0.431 0.347

179 80 -13 -103 166 80 -15 -105 164 73

6"351 2.577 1-536 1"070 0.798 0-599 0"439 0"312 0.215 0"146

-179 83 -8 -100 166 68 -32 -134 124 22

6"388 2.616 1"575 1-101 0"815 0.610 0.449 0-323 0"224 0.150

-178 87 4 -74 -149 140 71 5 -61 -125

6.403 2-608 1.525 1"000 0"673 0-448 0.290 0-183 0.112 0.067

-178 88 6 -72 -148 138 65 -7 -79 -150

6"345 2.574 1"540 1.084 0.831 0.664 0"538 0"430 0.334 0.247

-179 88 -8 -98 171 78 -18 -117 142 38

6"362 2"591 1-561 1-109 0.856 0.687 0"555 0.440 0"345 0.258

-179 86 0 -82 -161 123 51 -18 -85 -148

TABLE 12

Transmitted powers for case 4--fine mesh ( N N = 3 0 1 ) xl0 s n 1 2 3 4 5 6 7 8 9 10 Total

Myers

Scheme I v = 0.01 X = 1/3

S c h e m e II v = 0.0 X = 1/3

S c h e m e II v = 0.01 X = 1/3

Scheme I v = 0.01 X = 1"0

S c h e m e II v = 0.0 X = 1.0

6-32 1"07 0.38 0"19 0.11 0.07 0"05 0"04 0"03 0-02 8-44

6"48 1"07 0.38 0-18 0.10 0"06 0.03 0-02 0"01 0.00 8.32

6"55 1-10 0"40 0"19 0.11 0"06 0"03 0-02 0.01 0"00 8.47

6-58 1.09 0.37 0-16 0.07 0.03 0'01 0.01 0-00 0-00 8-34

6.47 1"06 0"38 0.19 0.11 0"07 0"05 0.03 0"02 0-01 8.39

6"50 1"08 0"39 0"20 0.12 0"08 0-05 0-03 0.02 0-01 8.48

ACOUSTIC CHOKING

99

viscous t e r m s t e n d to r e d u c e the shock s t r e n g t h a n d h e n c e t h e s h o c k fails to d i s s i p a t e as m u c h e n e r g y as it s h o u l d . T a b l e s 9 a n d 10 s h o w the F o u r i e r coefficients o f the a c o u s t i c velocity at x = l a n d t h e t r a n s m i t t e d p o w e r o b t a i n e d with s c h e m e s I a n d II o n a m e s h with 251 grid p o i n t s f o r case 4. T h e c o r r e s p o n d i n g results o b t a i n e d on a mesh with 301 n o d e s are p r e s e n t e d in T a b l e s 11 a n d 12. T a b l e s 11 a n d 12 also i n c l u d e the results o b t a i n e d with s c h e m e s I a n d I I on the 301 n o d e m e s h with X = 1.0. F i g u r e 9 shows the p l o t s o f a c o u s t i c v e l o c i t y given b y each o f t h e s c h e m e s o n this mesh. The a s y m p t o t i c s o l u t i o n is p l o t t e d with the n u m e r i c a l x l O -4 8 000 6.000

\

~.

.

(o)

4 000

"9

2"000

.~

0000

~~27

O0

28'00

2900

~50 O0

-2.000 "~ - 4 . 0 0 0 -6 000 -800C 8000 • .4 6 000

,•,

_--\

(b)

4 000

5.~o 2 0 0 0 >~ 0 0 0 0

~ -2"000 -4,ooo

L,

O0

28 O 0

ime

\

/

-6000 -8 0 0 0

Figure 9. Acoustic velocity at x=! for case 4. NN=301, g = l . 0 . (a) Scheme I, J,=0.01; (b) scheme II, v = 0.0. - - - , Myers; - - - - - , linear; , numerical solutions. TABLE 13

Fourier coej~cients for case 5--coarse mesh ( N N = 251, X = 1/3) xl03

Myers

Scheme I I,=0.1

Scheme I .v=0.05

Scheme II z, =0.01

Scheme II ~,=0.015

.

IbJ

0~,

IbJ

0~,

IbJ

0~,

Ib.I

0~,

Ib.I

0~,

1 2 3 4 5 6 7 8 9 10

0"998 0.495 0.327 0"244 0-195 0.162 0"139 0.122 0.108 0"097

172 48 -70 174 59 -56 -171 74 -41 -156

1"043 0.522 0.354 0.282 0-252 0"235 0.214 0"186 0.156 0.126

177 58 -53 -164 84 -32 -151 90 -31 -152

1-059 0"583 0"459 0-367 0-270 0-188 0.127 0.084 0"055 0.035

176 57 -55 -165 85 -27 -139 108 -6 -120

1"058 0"532 0"371 0.309 0.279 0.251 0.213 0.170 0"128 0.093

180 65 -44 -148 113 21 -68 -156 117 29

1"051 0"533 0"382 0.327 0"292 0-250 0"199 0.150 0"109 0.078

178 61 -49 -154 108 16 -75 -166 104 13

N. J. ~,VALKINGTON AND XV. EVERSMAN

100

TABLE 14

Transmitted powers for case 5--coarse mesh

(NN=251, x=l/3)

xlO s n

Myers

Scheme I u=O.O1

Scheme I v=O.05

Scheme II v=O.O1

S c h e m e 11 v=O.O15

1 2 3 4 5 6 7 8 9 10 Total

1.60 0.39 0" 17 0"10 0"06 0.04 0-03 0.02 0.02 0.02 2-57

1.75 0.44 0"20 0-13 0"10 0.09 0-07 0.06 0.04 0.03 2.94

1.80 0.55 0.34 0.22 0,12 0"06 0-03 0-01 0.01 0-00 3.12

1"80 0.46 0.22 0.15 0"13 0.10 0"07 0.05 0.03 0.01 3.03

1-77 0"46 0.23 0-17 0-14 0.10 0-06 0-04 0.02 0.01 3.01

TABLE 15

Fourier coefficients for case 5--fine mesh

x 10~

. 1 2 3 4 5 6 7 8 9 10

Myers

Ibol

o:

0.998 172 0.495 48 0-327 - 7 0 0,244 174 0.195 59 0.162 - 5 6 0.139 -171 0.122 74 0-108 -41 0-097 -156

(NN=301)

Scheme I u=O'Ol X=I/3

Scheme I v=O'05 X=I/3

Scheme II v=O-O1 x=t/3

Scheme II v=O.O15 x=l/3

Scheme I v=O-O1 x=l/3

Scheme II v=O-Ol X=I-0

Ib.I

Ib.I

Ib.I

lbol

Ibnl

Ibnl

o:

1.042 177 0.518 58 0.348 -53 0.267 -164 0-225 86 0-205 -25 0,194 -139 0.184 105 0,169 -14 0,150 -134

o:

1-048 177 0~541 57 0.405 - 5 6 0.348 -168 0.292 81 0,231 - 3 0 0.174 -142 0.127 105 0.092 -9 0.065 -123

0:,

1.059 180 0.528 65 0.358 - 4 2 0-282 -148 0.247 108 0.229 9 0,213 - 8 6 0.191 -178 0.162 93 0,132 4

0:,

1,052 179 0-526 62 6-360 - 4 8 0-291 -156 0-260 100 0.240 2 0-215 - 9 3 0,183 174 0,148 82 0.115 - 1 0

0:,

1.042 177 0-518 58 0.347 - 5 4 0,264 -165 0.216 84 0.185 - 2 8 0.164 -141 0,148 105 0.135 - 1 0 0,122 -129

07,

1,063 -179 0-528 68 0.354 -38 0,270 -142 0-222 115 0,191 15 0,169 -83 0,151 -177 0-134 92 0-115 4

TABLE 16

Transmitted powers for case 5--fine mesh ( N N = 3 0 1 ) X 107 n 1 2 3 4 5 6 7 8 9 10 Total

Myers

Scheme I u=O.O1 X=I/3

Scheme I v=O.05 X=1/3

Scheme II v=O.O1 x=l/3

S c h e m e II v=O.O15 x=l/3

Scheme I u=O.O1 x=l.0

Scheme.II v=O.Ol x=l.0

1 "60 0.39 0-17 0-10 0"06 0.04 0.03 0.02 0"02 0.02 2-57

1-74 0.43 0-20 0-12 0.08 0.07 0.06 0.05 0.05 0.04 2.91

1.76 0.47 0.26 0.19 0.14 0.09 0.05 0.03 0.01 0"01 3.02

1.80 0"45 0-21 0-13 0.10 0.08 0.07 0"06 0.04 0"03 3"01

1"78 0.05 0-21 0"14 0"11 0.09 0.07 0"05 0"04 0"02 2"98

1-74 0"43 0.19 0"11 0"08 0.06 0.04 0"04 0"03 0"02 2"79

1"82 0"45 0-20 0-12 0"08 0-06 0-05 0-04 0-03 0-02 2.88

ACOUSTIC CHOKING

101

solutions in Figure 9 and the solutions agree qualitatively. The results in Tables 9-12 show that the Fourier coefficients and the acoustic power a p p e a r to converge to the asymptotic coefficients and power; however the convergence is very slow. Tables 13 and 14 list the Fourier coefficients and power obtained with schemes I and II on a mesh of 251 nodes for case 5. Both schemes I and I I require a viscous term for this problem and the results obtained with two viscous coefficients for each scheme are presented. Scheme I is stable for relatively large viscous terms (for this C F L number); however scheme II does not enjoy this property. The results for case 5 obtained by these two schemes on a mesh of 301 nodes are presented in Tables 15 and 16. Comparison of the results for scheme I obtained with a viscous term o f ~, = 0.01 and those with v = 0.05 shows how an excess of viscosity can deteriorate the solution. Plots of acoustic velocity at x = l obtained with each of these two viscous terms are shown in Figures 10(a) and (b). It is difficult to tell by looking at the plots that the solution with v = 0.05 gives less

x10 "3 2500

(ol -%

2000

/\

"\

1500 ~>" 1000 o 0 500

o

~

0000

2

1~'27]00 i /

ooo

-0 500

,,~ - 1 0 0 0 - 1500

2000 1.5o0

O0

Time

\../

-2O0O -2000 x10-3

2500

\ /["

~'30

\

,-,\

(b) i

->" 1000 __o 0 5 0 0

0000 ~ -0500

,~

-

~

~ 3 0

O0

ooo

-1.500 -2000 -2500

Figure 10. Effect of the viscous term on the solution given by scheme I for case 5. NN=201, X = 1/3. (a) ~,=0"01; (b) v=O'05. accurate Fourier coefficients; however both give qualitative agreement with the asymptotic results. The acoustic velocity obtained with schemes I and II on the mesh of 301 grid points with X = 1.0 is plotted in Figure 11. T h e results for case 6 obtained by using schemes I and II on the 301 node mesh are given in Tables 17 and 18. Table 17 lists the Fourier coefficients of the acoustic velocity at x = l and Table 18 the transmitted power. The difference results obtained for this problem do not agree favorably with the asymptotic results also listed in Tables 17 and 18. The transmitted power given by scheme I is overestimated by a factor of about two and scheme II overestimates the transmitted power by a factor of four. The Fourier coefficients in Table 17 also exhibit similar errors. Table 4 gives a clue as to why such errors may be present. The asymptotic results predict that only 1.1% of the source energy emanates from the duct at x - - l. Scheme I then predicts that 2.4% of the source energy

102

N. J. %VALKINGTON A N D W. E V E R S M A N x10 -3

~

2.500 2000

15OO _~

1000

o

0 500

l

0o00

(o ) '

/-.

29 O0

~176 -0500

~ 3 0 O0

%

Time

-1 0 0 0 -1-500

,, /

-2000 -2 500 x I 0 -3 2500 2000

-~ 1-500

-~

/--.\

\

( b}

%

I ooo

._o 0 5 0 0

g 00oo (? - o 5 o o

2obc Time

- 1-000

\

-1.500

-2-000

/

,

i

\./

-2"500

Figure 11. Acoustic velocity at x = 1 for case 5. N N = 3 0 1 , - - - , Myers; - - - - - , linear; - - , numerical solution.

X =

1.0; v = 0 . 0 1 . (a) Scheme I; (b) scheme II.

is transmitted and scheme II gives 4.8% of the energy transmitted. The results for cases 3, 4, and 5, when expressed as a ratio of transmitted to source energy, give errors similar in magnitude to this. This suggests that the problems encountered with case 6 may be caused by the solution being similar in size to the errors in the scheme being used to produce it. Figure 12 shows the plots of acoustic velocity at x = 1 for case 6 obtained by using the 301 node mesh with X = 1.0. Comparison ofthe linear solution with the difference solutions shows that the transmitted wave is considerably smaller. However it is very apparent in T A B L E 17

Fourier coejJicients for case 6 ( N N = 3 0 1 )

•

4

Myers

n

[bJ

o.~

1 2 3 4 5 6 7 8 9 10

1"934 1"031 0'696 0"524 0"420 0"350 0"300 0"263 0"234 0"210

-44 -35 -19 -1 19 39 59 79 99 120

Scheme I v=O.Ol X=I/3

Scheme I v=O'05 x=l/3

Scheme II v=O'Ol x=l/3

Scheme II v=O'O15 x=l/3

Scheme I v=O'Ol x=l.O

Scheme II v=O-Ol X=I-O

Ib.I

[bJ

[b.I

Ib.I

[b.I

[b.I

o~,

2.743 -7 1"708 33 1.297 79 0"980 131 0-668 - 1 7 5 0"433 - 1 2 5 0"282 - 7 7 0"184 - 2 9 0"119 17 0'076 64

0~,

2"917 - 2 0 1"983 25 0-926 102 0"347 156 0-158 - 1 4 4 0-063 - 8 5 0"026 - 3 1 0"011 22 0-005 71 0.003 - 1 2 2

o~,

3.864 42 2"201 140 1.705 - 1 1 1 1-460 8 1.185 135 0"856 - 9 3 0.574 36 0-381 164 0"255 - 6 7 0-169 63

o,~

3.874 41 2-275 137 1"803 - 1 1 4 1"490 3 1-086 137 0"708 - 9 5 0-454 31 0"294 157 0"190 - 7 5 0"119 54

o:

2"614 -8 1"455 34 1-109 90 0-935 139 0-787 - 1 7 3 0"635 - 1 2 4 0"483 - 7 5 0.347 - 2 8 0.240 18 0-160 63

o~,

3"766 41 2-027 137 1"436 - 1 1 7 1"173 -8 1.027 107 0-901 - 1 3 2 0"748 -5 0"567 126 0"392-102 0"252 30

ACOUSTIC

103

CHOKING

TABLE 18

Transmitted powers for case 6 ( N N = 301)

gl

Myers

Scheme I v=0.01 X= 1/3

1 2 3 4 5 6 7 8 9 10 Total

0"59 0.17 0.08 0.04 0"03 0.02 0"01 0.01 0,01 0"01 1.03

1-20 0-47 0.27 0"15 0-07 0"03 0"01 0.01 0-00 0.00 2.21

xlO 8

Scheme I v=0.05 X=1/3

Scheme II v=0.01 X= 1/3

Scheme II v=0,015 X= 1/3

Scheme I v=0.01 x=l.0

Scheme II v=0.01 X=I.0

1.34 0"63 0"13 0"02 0"00 0"00 0.00 0-00 0"00 0"00 2.15

2-39 0-77 0"46 0.34 0.22 0" 12 0.05 0.02 0-01 0.00 4.40

2"40 0-83 0.52 0.35 0.19 0.08 0"03 0"01 0"01 0.00 4-42

1-09 0.34 0.20 0.14 0. I0 0"06 0-04 0"02 0.01 0.00 2.00

2.27 0.66 0.33 0"22 0" 17 0" 13 0-09 0.05 0"02 0"01 3.95

x10-3 2500 2000

./-\,

i ol

\

1.500 ~" 1 0 0 0 o 0500 ,.J 0 0 0 0 ~ -0500

/

~ -looo - 150O -2000 - 2500 x l 0-3 2500

-~ 0

r, me

\

/

i

!

\ i (b)

Z000 1500 1000 05o0

.e 0 0 0 0 '~ - 0 5 0 0

<~- ~ ooo - 1500

-ZOO0 -2500

"V---b--'

/

/

2 oo

,ime \

,_1"

F i g u r e 12. A c o u s t i c v e l o c i t y at x = l f o r case 6. N N = 3 0 1 , - - -, Myers; .... , linear; , numerical solution.

',,

/

/

\.j X = 1.0, v = 0 . 0 1 .

(a) S c h e m e I; (b) s c h e m e II.

Figure 12 that the shock location given by the differencing scheme does not coincide with the location predicted by the asymptotic results. If the shock wave is not calculated correctly the power calculations should be in error since it is the shock that dissipates the energy. One error always present in the differencing schemes is the drift of the mean flow from the steady mean flow established prior to the introduction of the acoustic signal. For cases 1 to 5 the drift in the mean flow velocity was at least two orders of magnitude smaller than the fundamental acoustic harmonic amplitude of the source. However, for

104

N. J. "~VALKINGTON A N D x,V. E V E R S M A N

case 6 the transmitted signal is also two orders o f magnitude smaller than the source amplitude. The asymptotic theory does admit the possibility of acoustic streaming, which is a change in the mean flow caused by the acoustic signal; however for all of the problems solved here the streaming at the duct terminations is negligible. When the throat Mach n u m b e r is close to unity a small change in the steady Mach n u m b e r at the duct terminations will cause a large change in the throat Mach n u m b e r (see the p a p e r by Myers and Calegari [5]). A numerical experiment was performed to see if adjusting the initial conditions so that the mean throat Mach number in the last period was correct to four decimal places would improve the solution. The results obtained with this procedure by using scheme I are presented in Table 19. The Mach numbers at the duct inlet, outlet and throat are also listed in Table 19 prior to introducing the acoustic signal and in the last period. The results obtained by using this procedure are no better than t h o s e obtained previously. TABLE 19

Results of numerical experiments for case 6 Experiment 1

Experiment 2

n

]b~lx 10'

0,~

Power

(b,l• 104

0:

Power

1 2 3 4 5 6 7 8 9 10

2.916 1.589 1.185 1-001 0.863 0.716 0.563 0.418 0.297 0.202

6 66 132 -164 -102 -40 22 83 142 -159

1-36 0-40 0-22 0.16 0.12 0.08 0-05 0.03 0.01 0.01

2.149 1-258 0.973 0-800 0-644 0.487 0-344 0.233 0-153 0.097

-31 -10 13 36 61 87 111 133 154 173

0.74 0.25 0-15 0.10 0.07 0.04 0-02 0.01 0-00 0.00

--

2.45

--

1.38

Total powerx I0 s Throat Mach numbert Throat Mach number:l:

-0.978358 -0.980070

-0.982871 -0.983223

Exit Mach humbert Exit Mach number~:

-0.508658 -0.508575

-0-508765 -0-508677

t Steady conditions prior to acoustic signal. ~t Mean conditions in final period.

Table 19 shows that the steady throat Mach number prior to introduction of the acoustic source is -0.978. This throat Mach number being further from - 1 . 0 than - 0 . 9 8 will tend to decrease the non-linear trends. Table 19 also lists the results obtained with scheme I when the initial conditions were adjusted to give a steady mean flow throat Mach number of -0.983 prior to the introduction o f the acoustic signal. The tabulated powers for this example show that this slight increase in Mach n u m b e r has reduced the transmitted power significantly. Figure 13 shows the plot of acoustic velocity at x = ! for these last two experiments and it can be seen that the latter gives a shock position similar to that predicted by the asymptotic results. In summary, the results for this case 6, which has a throat Mach number very close to unity, are very sensitive to the mean flow profile since this profile will influence the shock strength. Accordingly small numerical errors may lead to poor results. Such errors will

ACOUSTIC CHOKING x10-3 2 500 2 000 1.500 o o

-~

,,.'

- 1500

-2000 -2 500 xi0-3 2 500 2000

/

-

/

',.

(a}

0500 0 000

8

-0 500 -I 000 "-1.500 -2 0 0 0 -2.500

-'

\

ooo

28 O0 Time

\,

\

,../

/

/

.Y

/ "\

I "500 1000

= <

\

1 ooo

-~ooo

u

~-~

o 500

.-'=' 0 0 0 0 ~ ,~ - 0 500

~, ~

/

105

i

(b)

i

27oo// ""'~"/-~-'

,j/'

/

%

2900 3000

_

,. . . . . . .

Z8 O0

,me\

i

/

\,\./

Figure 13. Results of the numerical experimentswith case 6 Obtained by using scheme I. N N = 301, X= 1-0, v = 0-01. (a) Effectof decreasing the throat Mach number; (b) effect of increasing the throat Mach number. not only modify the mean flow profile but may be similar in magnitude to the acoustic perturbations. 5. CONCLUSIONS In aircraft turbofan inlets, fan generated noise is observed experimentally to be significantly attenuated at high subsonic inlet Mach numbers. This phenomenon cannot be predicted by linear theory. In order to study the physical process by which this may occur a numerical algorithm has been developed to solve a related non-linear problem in one dimensional gas dynamics. The non-linear solution predicts wave steepening and shock waves that do not appear with the linear theory. Such discrepancies may explain the high attenuations obtained experimentally. The acoustic solution is obtained by solving the transient gas dynamic equations. The acoustic quantities are evaluated by subtracting the steady state solution of the gas dynamic equations from the transient solution. Approximate solutions are obtained by using several finite difference schemes. All of the differencing schemes have been used to solve multidimensional problems in gas dynamics. The boundary conditions required to model an acoustic source and an anechoic termination developed previ6usly were found to be satisfactory for the shocked problem. Comparison of the results obtained with the different schemes indicates that they all give approximations to the correct solution with no one o f them giving significantly "better" results. The numerical algorithms duplicate the wave steepening and shock formation phenomena predicted by analytical methods. Decomposition o f the acoustic velocity into its Fourier components enables the numerical solutions to be compared with asymptotic results. The Fourier components of the acoustic velocity ~it the duct exit appear to converge to those given by asymptotic results as the grids are refined, however convergence is very

106

N. J. W A L K I N G T O N

AND

%V. E V E R S M A N

slow. Exact agreement between the numerical and asymptotic results is not expected since both methods o f solution are approximate. In order to model shock waves it is necessary to modify the basic differencing schemes. The most successful modification found was the addition of an artificial, fourth order viscous term. Numerical experiments suggest that this term changes continuous solutions negligibly. The viscous term improved the accuracy of the discontinuous solutions. Numerical experimentation is required to estimate a reasonable viscous coefficient since too large a coefficient will degrade the accuracy of the solution. Plots o f the acoustic velocity show that the shock locations predicted by the numerical schemes coincide with those predicted by the analytical solutions. The reduction in the ratio of transmitted to incident p o w e r that occurs with the shocked solutions indicates that the shock waves dissipate acoustic energy. This otters an explanation of the acoustic choking observed experimentally. The results indicated that more power is dissipated as the Mach number, sound amplitude and frequency are increased. These observations are in agreement with those observed experimentally and those predicted by the method of matched asymptotic expansions. Reflection of acoustic energy, back to the source, could explain the acoustic choking; however the non-linear solutions do not produce more reflection than the linear solutions. The one dimensional model used for this study can accommodate only plane waves. The multi-dimensional problem will have modal interactions which may significantly attect the non-linear solution. The one dimensional results suggest that the a p p r o a c h taken here will yield solutions to the multi-dimensional problem; however such a code will require prohibitively large processing times with the current generation of computers.

ACKNOWLEDGMENT The authors appreciate the assistance of Dr M. K. Myers who generously provided the asymptotic results for the Tsien-Crocco duct. This work has been supported by NASA Lewis Research Center grant N S G 3231. I M S L routine F F T R C was used to decompose the acoustic signals into their Fourier components. REFERENCES

9

9

1. N.J. WALKINGTON and W. EVERSMAN 1983 Journal of Sound and Vibration 90, 509-526. A numerical model of acoustic choking, part I: Shock free solutions. 2. J. VON NEUMANN and R. D. RICHTMYER 1950 Journal of Applied Physics 21, 232-237. A method for the numerical calculation of hydrodynamic shocks9 3. D. T. BLACKSTOCK 1961 Journal of the Acoustical Society of America 34, 9-30. Propagation of plane sound waves of finite amplitude in nondissipative fluids. 4. D.T. BLACKSTOCK 1966 Journal of the Acoustical Society of America 39, 1019-1026. Connection between the Fay and Fubini solutions for plane sound waves of finite amplitude. 5. M. K. MYERS and A. J. CALLEGARI 1977 Journal of Sound and Vibration 51, 517-531. On the singular behavior of linear acoustic theory in near sonic duct flows. 6. A. J. CALLEGARI and M. K. MYERS 1977 American Institute of Aeronautics and Astronautics Paper 77-1296. Nonlinear eiiects on sound in nearly sonic duct flows. 7. A. J. CALLEGARI and M. K. MYERS 1979 American Institute of Aeronautics and Astronautics Paper 79-0623. Sound in ducts containing nearly choked flows. 8. M. K. MYERS and A. J. CALLEGARI 1980 American Institute of Aeronautics and Astronautics Paper 80-1019. Acoustic shocks in nearly choked flows. 9. M. K. MYERS 1981 American Institute of Aeronautics and Astronautics Paper 81-2012. Shock development in sound transmitted through a nearly choked flow.

ACOUSTIC CHOKING

107

10. G. A. SOD 1978 Journal of Computational Physics 27, 1-31. A survey of several finite difference methods for systems of nonlinear hyperbolic conservation laws. 11. M. VINOKLIR 1974 Journal of Computational Physics 14, 105-125. Conservation equations o f gas dynamics in curvilinear coordinate systems. 12. G . W . HEDSTROM 1979 Journal of Computational Physics 30, 222-237. Non-reflecting boundary conditions for nonlinear hyperbolic systems. 13. A. HARTEN and G. ZWAS 1972 Journal of Computational Physics 9, 568-583. Self adjusting hybrid schemes for shock computations. 14. J. P. BORIS and D. L. BOOK 1973 Journal of Computational Physics 11, 38-69. Flux-corrected transport. I. SHASTA, A fluid transport algorithm that works. 15. D. L. BOOK, J. P. BORIS and K. HAIN 1975 Journal of Computational Physics 18, 248-283. Flux-corrected transport II: Generalizations of the method. 16. J. P. BORIS and D. L. BOOK 1976 Journal of Computational Physics 20, 397-431. Flux-corrected transport III: Minimal-error FCT algorithms. 17. R. W. MACCORMACK 1969 American Institute of Aeronautics and Astronautics Paper 69-354. The effect of viscosity in hypervelocity impact cratering. 18. G. MORE-rrl 1979 Computers and Fluids 7, 191-205. The Lambda scheme. 19. G. MORETTI and L. ZANNETI'I 1982 American Institute of Aeronautics and Astronautics Paper 82-0168. A new improved computational.technique for two dimensional unsteady compressible flows. 20. B. VAN LEER 1979 Journal of Computational Physics32,106-136. Towards the ultimate difference scheme. V. A second order sequel to Godunov's method. 21. N. J. WALKINGTON 1983 Ph.D. Thesis, University ofMissouri-Rolla. A numerical model for subsonic acoustic choking. 22. J. L. STEGER and R. F. WARMING 1981 Journal of Computational Pilysics 40, 263-295. Flux vector splitting of the inviscid gas dynamics equations with applications to finite difference methods. 23. G. D. VAN ALBADA, B. VAN LEER and W. W. ROBERTS 1982 Astronomy and Astrophysics 108, 76-84. A comparative study of computational methods in cosmic gas dynamics.