-16-

2
cl R
ops tf ap
5 ar (p - Pp) + “oe: (2.53)

The quantity |1/ol is, of course, given to us by Figure 1. The skin
friction coefficient could in principle be obtained from the actual
computation at any time level. However, within our approximation one may give
an estimate for Ce taken, for example, from boundary layer theory. Thus,
for laminar flow and subsonic Mach numbers we may use the Blasius value of

C, 2 2/3VR, . More on this in the section describing the numerical work
which tested the efficacy of the boundary condition (2.53).

We close this section by noting that the same methodology used here to
derive outflow boundary conditions for the flat-plate problem can be used to
extract boundary condition on the downstream side of other configurations such
as the conditions in the wake of an airfoil, flow in a pipe, and so on. Also,
other quantities have perturbations which satisfy (2.53), and one might try to
apply the outflow condition to them as well. Another potential application is
for unsteady flows where the dispersion relation for (2.13) - (2.15) would be

carried out around some central wavenumber By rather than around 8 = 0.

3. NUMERICAL RESULTS

In this section we present some numerical results for the boundary
condition described above. These results are preliminary; more complete
results will be presented elsewhere. Using the non-dimensionalization
described above and the Blasius expression for Cg, the boundary condition

(2.53) can be written as


-10-

(0) (0). _

dug 61°? + 65°? = 0 (2.29)

(0) (0) _ lO n
iw Gy“ + G3 Gy =0 (2.30)

(0) 1 (0). _ y-1 (0)\- (On .
tug oS? + of - - (vs: ) - 65" = 0 (2.31)
YM,
and
inp of) + of) = - to, of) - ic6) - ing (2.32)
(1) 4 6D _ gn LL n\g(0) _ i ,(0)

jug Gy" + Gy eG, = (6 + i)°% 2 Gy (2.33)

ivy of) + ar 6) - a (<) - eon = - io, + n)6,°, (2.34)

co

Similar sets of equations can be derived for higher powers in 8. They
all have the same left-hand side in the relevant variables, with right-hand
sides serving as source terms made up of all previous (presumably known)
approximations. It is seen that YQ» the zeroth-order approximation to the
decay rate is determined from the "8 = 0 problem," equations (2.29) - (2.31),
without reference to w= (du /d8),_9, the phase velocity of the perturba-
tions. We, therefore, first concentrate our attention on this "B=0
problem." By some simple manipulations one can extract from the system (2.29)
- (2.31) a single fourth-order differential equation for of, For the sake
of neater notation set of) = 85> (j = O,l,ee*). With this change of
notation, the fourth-order differential equation for & = of) is found to

be


-20-

It is apparent from Table II that the values for U are very insensitive
to the position of the outflow boundary’ and to the choice of (1.1) or (3.1).
The values obtained from (3.1) are slightly less accurate, but the differences
are very small. At grid point 15 the relative error is only about 1.5%.
These errors decrease for grid points away from the wall.

The results in Table III show that there are significant boundary errors
in V for both boundary conditions; however, these values are very small
compared to the free stream values. At grid point 13 the relative error in V
is about 8.5% for (3.1) and about 10.9% for (1.1). It would require a more
complicated problem to assess the meaning of these differences.

The boundary condition (3.1) has also been used to compute the steady
flow field over a curved surface. A particular example is the flow over a
curved hump with the curvature chosen so as to accelerate the flow from
M, = +7 to M, = 0.76. The height of the hump was .04 feet and the
curvature extended over a range of 1.2 ft. The surface then became flat and
the computational domain extended a further 1.2 ft downstream. The unit
Reynolds number was 3.0 x 10° while the boundary layer thickness at inflow was
+0095 ft. The Reynolds number based on the boundary layer height at inflow
was 896.00. Using (1.1), and several different values of a, convergence to
the tolerance described above could not be obtained. Using (3.1) convergence
in 52,760 iterations was obtained.

Finally, we point out that the methodology described above can be applied
to schemes that are inconsistent in time. As one example, we can consider

systems of the form

io]
|
+

w]e
ole]
aka
* ra]

ele
213;
at
Oo


1. INTRODUCTION

Since the early 1970°s there has been growing attention paid to the
numerical solution of the full time-dependent, viscous, compressible Navier-
Stokes equations of fluid mechanics; see inter~alia [1], [4]. A vexing and
still open issue is the formulation of "downstream" or "outflow" boundary
conditions. One common strategy for flows evolving to steady state has been
to assume that the particular geometry considered allows the imposition of the
infinity-downstream steady state conditions at the relevant boundary of the
computational domain. For example, for flow in a pipe, the downstream
condition might be taken from the fully-developed, steady-state parabolic
profile. For a semi-infinite flat plate at zero angle of attack, it is
usually assumed that "sufficiently far" downstream, the free stream conditions
predicted by boundary layer theory are sufficient for stating the
corresponding boundary conditions.

‘Numerical experiments (see e.g., [5-6]) with time-consistent codes have
shown that taking p = p,, for the downstream boundary conditions leads to
very slow convergence to steady state. This phenomenon is aggravated with
increased Reynolds numbers, because the spatial resolution of the boundary
layer then leads, via the von Neumann stability analysis, to very small time
steps. It can be seen from the results of Rudy and Strikwerda [5] that the
retardation of the convergence to steady state is caused mainly by the slow
decay of relatively long wavelength disturbances. They tried to account for

this by proposing a boundary condition on the pressure of the form

BP - pc M2 + a(p - p,) = 0 a1)


-17-

op. - ap
38 A(p - p,) + B 3 (3.1)

where

We implemented (3.1) in a Navier-Stokes program using a version of the
MacCormack scheme that is second-order accurate in time and fourth-order
accurate in space [7]. The numerical scheme was time consistent and did not
use any acceleration technique such as a local time step.

For our numerical experiment we used an M, = 0.4 boundary layer. The
unit Reynolds number (per foot) was 3.0 x 10°. The inflow was taken at 1.0
feet from the leading edge. We performed calculations using two different
locations of the outflow boundary; x = 2.067 feet and 3.0 feet from the
leading edge respectively. In both cases we used a grid of 31 points in the y
direction (exponentially stretched). In the x direction we used 17 and 31
grid points respectively for a mesh size of .067 feet. The inflow and initial
data were obtained from a computer program which solved the boundary layer
equations. Convergence to steady state was assumed when the maximum over the
grid of 22 was less than 107°,

The boundary condition (3.1) was compared with the condition (1.1).
Equation (3.1) has the property that the steady state pressure is not forced
to be equal to p. Thus, there may be some discrepancies in the steady state
solution. The numerical experiments were designed to assess the behavior of

(3.1) with regards to:


[8]

(9]

-23-

Harris, J. E. and D. K. Blanchard, Computer Program for Solving Laminar,
Transitional or Turbulent Compressible Boundary Layer Equations for Two-

Dimensional and Axisymmetric Flow, NASA TM 83207 (1985).

Turkel, E., Algorithms for the Euler and Navier-Stokes Equations for

Supercomputers, NASA CR _ 172543, Progress and Supercomputing in

Computational Fluid Dynamics, E. M. Murman and S. S. Abarbanel, (eds.),

Birkhauser Boston, Inc., 1985, pp. 155-172.


~19-

case of the outflow boundary at 3.0 feet. In this case the outflow pressure
obtained from using (3.2) differed from P,, by less than .04%. In Table I we
indicate the value of U and V_ (horizontal and vertical velocities non-
dimensionalized by U,) for the steady states obtained from (3.1) and (1.1).
The data are presented at different x locations indexed by the grid point.
(The outflow was at grid point 31.) In all cases the velocities are shown for
the first grid point away from the wall as this was where the maximum relative
difference occurred. The relative differences decreased significantly for
grid points away from the wall where the velocities increased in magnitude.

Using (1.1), we observed a slight oscillation in V_ near the outflow.
This may be due to the fact that the outflow pressure is fixed at the free
stream pressure while the flow appears to undergo a very small acceleration.
For the boundary condition (3.1), V decreases rapidly close to the outflow.
Thus, V near the outflow appears to be incorrect with both boundary
conditions. Away from the outflow it is clear that V is the quantity that
is most sensitive to the boundary condition although the differences are small
and may also be affected by truncation error. The horizontal velocities are
very close. At grid point 27 the two values of U differ by about 1.7%.

In Table II we list the values of U obtained from using (1.1) and (3.1)
with the outflow boundary at 2.067 feet from the leading edge (grid point
17). These are compared with the values obtained when the outflow boundary
was at 3.0 feet from the leading edge. We assume these values are accurate
for the problem of the plate extending to infinity. Ideally, one would want
the solution to be independent of the position of the outflow boundary. The
values in these tables are the first grid point away from the wall where the

relative differences are greatest.


-14-

where in (2.40) - (2.46) for given o ande we already know, albeit
numerically, gq and Wye In deriving the boundary conditions (2.45), (2.46)
one has to expand the expression (2.19) for r in powers of 8. One then
finds that the expansion is valid only for Bly small compared to unity.
(Exactly how small cannot be answered at this point without finding the radius
of convergence for our expansion.) In physical terms, for laminar flow where
é~ vx, this constraint means that our analysis is restricted to modes whose
wavelength is longer than the computational domain on the plate.

A solvability condition for the 8)~problem is found by multiplying (2.40)
by the left null vector of the left-hand side operator of (2.40).

Symbolically, we may write

<£,R,>

1.77 ED (2.47)

where f£ is the left null vector to the system, including the boundary
conditions. Numerical evaluation of (2.47), for the same range of the
Reynolds and Mach numbers mentioned before, leads to the result that wy is

almost a universal constant in the present problem, namely

w) = — 2295 + .002. (2.48)

Now that we have the decay rate w and the phase velocity w) for long

0

waves, we can derive outflow boundary conditions.

Consider, e.g., the pressure perturbation:


Thus, p’ is the pressure perturbation, etc. Next, we linearize the Navier-
Stokes equations in the boundary layer region (0 < y < 6) by substituting
(2.7) into (2.1) and retaining only first-order quantities in the primed

variables. The resulting set of equations may be written as follows:

op7 op7

ap TAY Gt Pol5x + by ) = 0, (continuity) (2.8)

au’ au - 1 apt _

Pra + ay x + av? + box Yo 7 > (x-momentum) (2.9)
0 ay

Ov ggy yl aly aw (y tum) (2.10)

at Y ox 65 ay % ay > y-momentum) , .

and instead of the linearized energy equation we have the linearized constant

enthalpy relationship:

Lo)

=v 0 - -
=I Po = yoT ae Pg Ug ue (2.11)

In the above equation Yo is the kinematic viscosity at steady state
conditions. We have also used the thin layer approximation in (2.9) and
(2.10). We emphasize again that the perturbation equations (2.8) - (2.11) are
valid only for O<y< 6. Outside this viscous region we shall derive
another set of equations which then will have to be matched to the above
equations. First, however, we set the modal analysis. Because (2.11)
provides an algebraic relationship between p” and p” and =u’, it is

necessary to consider only a three-vector, (p”, u”, v~) for which we set the

modal ansatz:


2, Government Accession No.

” Report No. NASA CR-178068
ICASE Report No. 86-9

3. Recipient's Catalog No.

4. Title and Subtitle

NON-REFLECTING BOUNDARY CONDITIONS FOR THE COMPRESSIBL!
NAVIER-STOKES EQUATIONS

5. Report Date
March 1986

7. Author(s)
Saul Abarbanel, Alvin Bayliss, and Liviu Lustman

6. Performing Organization Code

8. Performing Organization Report No.

86-9

9. Performing Organization Name and Address
Institute for Computer Applications in Science

and Engineering
Mail Stop 132C, NASA Langley Research Center
Hampton, VA 23665-5225

10. Work Unit No.

11. Contract or Grant

No.
WAS1-17070, NAS1-18107

12, Sponsoring Agency Name and Address

National Aeronautics and Space Administration
Washington, D.C. 20546

13. Type of Report and Period Covered

Contractor Report

14. Sponsoring Agency Code

505=-31-83-01

15. Supplementary Notes

Langley Technical Monitor:
J. C. South

Science &
Final Report

Invited Lecture, 7th Intl. Conf.
on Computing Methods in Applied

Engineering, 12/9-13/85,

Versailles, France (INRIA)

16. Abstract

A small perturbation analysis, in the long wavelength regime, is used to
obtain the downstream boundary condition for the pressure for the flow over a

flat plate.

results for high Reynolds number laminar flows

The methodology is extendable to other geometries.
show great improvement in

Numerical

convergence rate to steady state as well as in the quality of the results.

17. Key Words (Suggested by Authors(s))

boundary conditions
Navier-Stokes equations

18. Distribution Statement

34 - Fluid Mechanics & Heat Transfer
64 - Numerical Analysis

Unclassified ~ unlimited

19. Security Classif.(of this report)
Unclassified

20. Security Classif.(of this page)
Unclassified

21. No. of Pages] 22. Price
27 A03

For sale by the National Technical Information Service, Springfield, Virginia 22161

NASA-Langley, 1986



2. THEORETICAL DEVELOPMENT

The geometry under consideration is shown in the sketch below.

Computational domain
Inflow boundary tow bound
NQpce tec cease scott ow boundary

1 edge" of

boundary layer

i

y —_—_—_—

L..

1
{
'
|
4
1

The two-dimensional compressible Navier-Stokes equations governing the flow

may be written in conservation form as

aU | OF , ac
oh 4 224 Se .
xt +3 By (0) (2.1)
where
p pu
2
au: pu Fe pu + Tex
pv , puv + T
yx at
E (E+ To)" + Ty V7 k ox
(2.2)
pu
uv + T
Ge ° xy
pv2 +T .
y ar

(E + Ty) + Tay u-k ay


where p is the density, c = (yp/p) Mp is the speed of sound, and a is a
free parameter which is "tuned" so as to optimize the convergence rate to
steady state. Note, that in any case at steady state, the downstream boundary
pressure has the free-stream value. Their treatment greatly accelerated the
convergence rate; its derivation, however, is basically "hyperbolic."

In this paper we propose to derive downstream boundary conditions using a
different methodology. The basic idea is to perturb the Navier-Stokes
equations around some approximation (usually a fairly crude one) to the steady
state solution in the downstream region. The linearized partial differential
equations for the perturbed quantities are then attacked using a modal form
for their solution. The resulting ordinary differential equations present an
eigenvalue problem which is solved assuming very long wavelengths. This
eigenvalue problem yields the decay rate for long waves and their phase
velocity. With these two quantities, as functions of the Mach number and the
Reynolds number, one can use the modal form of the perturbed solution to
obtain the desired boundary condition.

In section 2 we give the derivation of the eigenvalue problem and the
boundary condition for the particular case of the flow past a two-dimensional
flat plate.

Numerical results are presented in section 3 for the case of subsonic
Mach number, M, = 0.4. The corresponding Reynolds numbers per foot
is 3x 105, The numerical experiments are carried out for semi-infinite
plates with the outflow boundary taken at various locations downstream of the

leading edge.


-15-

pe et(yt + bx)

G(n) (2.49)
where G is a linear combination of G,; and Gy. Since the numerical
solution of the finite difference approximation to the Navier-Stokes equations
is often done in terms of dimensionless temporal and spatial variables, we

define dimensionless time and x coordinates by

(2.50)

o
iy
om
4

alts

where the reference length L is the distance from the leading edge to the
outflow. Using y = aw, b = af/U,, and (2.50), we can cast (2.49) in the form
R, ¢
co

i> £ (we + BE)
po~e G(n) (2.51)

where Ce is the skin-friction coefficient,

1 =
Te = . (2.52)
Using the long-wave dispersion relation

0 =o, + Buy

0

and p° = p-p,,, we find that


We assume, first of all, that near the outflow boundary the steady state
is given by a constant pressure Pg» constant density Por a vanishing
vertical velocity vg = 0, and that the velocity in the x-direction is given
by
(6<y<-)

Ug = (2.5)
ay (0 <y< 6)

where 6 is the nominal edge of the boundary layer and a= t/y is the
velocity gradient at the wall (y = 0). The shearing stress, t, is, of course,
a function of x but will not be treated as such in the development.

A second simplification is to assume (for the mean steady flow) that the

total enthalpy, Hg, is constant, i.e.:
2
vara tz pty ) = H, = const. (2.6)

Equation (2.6) is equivalent to the assumption that we are perturbing around a
steady flow of Prandtl number unity.

The perturbed field quantities are defined as departures from the assumed
steady state described above:

uw=u-u

0

Po =P ~ Pg
(2.7)

P7 =P -


~21-

where the flux functions F and G are the same as in (2.1) and E(U) is a
preconditioning matrix (see [9] for appropriate examples). The matrix E(U)
can be included in the linearization with only slight changes in the algebra.
ay rg

(Note that is zero so that the only modification is that oe

3t _ is
replaced by E(U,) au")
O° ot °°

A simpler example is the use of a uniform grid in the viscous region and

a different (larger) uniform grid in the inviscid region with the local time
step appropriate for stability used in each region. It is easy to see that on
the differential equation level this can be modelled by a matrix E which is
I in the viscous region and of the form tt where B is some large number
depending on boundary layer thickness. This only affects the inviscid
matching conditions ((2.22) - (2.24)). We have found that the parameters A

0
and w, are insensitive to B, if B= ote).

4. CONCLUSION

We have introduced a methodology which enables us to derive an outflow
boundary condition for a range of steady flows. Preliminary results indicate
that the resulting boundary condition may be effective in accelerating the
computation of steady subsonic flows with little loss of accuracy and no free
parameters to be determined. The method is also applicable to incompressible

and supersonic flows. Extensions to unsteady flows are also possible.


NASA ~CK-V18 COS
NASA Contractor Report 178068
ICASE REPORT NO. 86-93

ICASE

NON-REFLECTING BOUNDARY CONDITIONS FOR THE
COMPRESSIBLE NAVIER-STOKES EQUATIONS

NASA-CR-178068
19860013422

Saul Abarbanel
Alvin Bayliss

Liviu Lustman

Contract Nos. NASI-17070 and NAS1~18107
March 1986

INSTITUTE FOR COMPUTER APPLICATIONS IN SCIENCE AND ENGINEERING
NASA Langley Research Center, Hampton, Virginia 23665

Operated by the Universities Space Research Association

NASA LEBBARY GORY

National Aeronautics and APR 22 1986
Space Administration
Langley Research Center LANGLEY Resta | CENTER

LIBRARY, MASA

Hampton, Virginia 23665
HAMPTON, VIRGINIA


idw + B)G, + i8G, + 63 =0 (2.16)
B ip
iw+5e,+"c¢ =0 (2.17)
yr re 1
iv +e)¢, ++ eo -Yle =0 (2.18)
3 201 Y 2 ° °
YM,
rm

The solvability conditions for the existence of solutions of the form e
for the above homogeneous set of differential equations is found, after some

manipulation, to be

P= 8? - yt w + So + 8). (2.19)

The set of equations (2.13) - (2.15) constitutes a fifth-order system and
requires five boundary conditions. Two of them are supplied by the no~slip

condition at the plate:

6y(0) =

fl
°

(2.20)

G3(0) = 0. (2.21)

The other three are supplied by requiring the correct "impedance-matching" at

n = 1 between the system (2.13) - (2.15) and the set (2.16) - (2.18), namely:

671) = rG,(1) (2.22)

63(1) = r6,(1) (2.23)

63()) = r6,(1), (2.24)


7 Fity)
uv f= el(yttbx) Fo(y) . (2.12)
ve FAM)

Next, we normalize and nondimensionalize the various quantities as follows:

3
Q
r

ay/U, 17 Fyleas G, = F,/U,, G, = F,/U,,

wo = v/a, B = U,b/a,

o
0
g
<
2
a
sn
4
e
~
P=

where Ry is the Reynolds number based on the boundary layer thickness.
In terms of the dimensionless variables and parameters, after

substituting (2.12), equations (2.8) - (2.10) become, for 0 <n < 1:

i(u + Bn)G, + 486, + 6; = 0 (2.13)
8B iB = ect
idw + y n)G, + =z 6 + G = eG) (2.14)
YM,
i@w + ene, ++. o7 - Yb nes - 1G, = con (2.15)
3 ye 1 oy Me 2 = &63+ :

0

In equation (2.13) - (2.15) the primes designate differentiation with respect
to the coordinate n; e.g., 63 = dG,/dn, G3 = 2 G,/dn?.

Next, we consider the inviscid region, n > 1 where Up = U,, and the
flow is described well by the Euler equations (ce = 0). Repeating the above

procedure the acoustic perturbation equations corresponding to (2.13) - (2.15)

are derived:


‘ni wil i iu

3.1176 013!

DO NOT REMOVE SLIP FROM MATERIAL

Delete your name from this slip when returning material

to the library.

NAME DATE

MS

NASA Langley (Rev. Dec. 1991) RIAD N-75





and
du du, av
Tee Po 2u aE 7 ME * By)
= av du. ov
Tyy 7 P > 2H ay Gx * By)
= = - (2% 4 2%
Ty 7 Tyx = “Gy + Ix)

The various parameters and dependent variables P, uy, v, wy A, E, k, T, and
P are, respectively, the density, the x and y components of velocity, the
shear and second coefficients of viscosity, the total energy per unit volume,
the coefficient of heat conduction, the temperature, and pressure, An
equation of state p= p(e,p) relates the pressure to the density and the

specific internal energy e,
e = (E/o) ~ 4 (u? + v4), (2.3)
For example, for ideal gases the equation of state is
P= (y-L)ele +5 (u? + v9] (2.4)

where y is the ratio of specific heats at constant pressure and volume,
respectively.

The numerical algorithm used in the next section is based on solving a
set of difference equations which approximate the full Navier-Stokes equations
(2.1) - (2.4). However, in considering the perturbation model we make several

rather crude simplifications.


-18-

a) Convergence rate
b) Robustness with respect to changes in the coefficients

c) Quality of the steady state.

We summarize our results below.

Convergence to steady state was significantly accelerated by the use of
(3.1). The boundary condition (1.1) requires specification of a parameter a.
The convergence rate is sensitive to the choice of a. In [5] an optimal a
for a class of problems is presented. Because of the length of the computer
runs, we were unable to determine the optimal a for the present problen.
For the case of the outflow boundary at x = 3.0 feet from the leading edge a
steady state was obtained using (3.1) in 30880 time steps. Using (1.1) and
a= 1.45 the steady state required 97020 time steps. Using a = .6 which
would be obtained from the formula in [6] based only on the free stream
velocities and sound speed significantly degraded the convergence rate. For
the case of the boundary at x = 2.067 feet from the leading edge, the use of
(3.1) required 20960 time steps for convergence. The same run with a = 2.71
required 52080 time steps.

Robustness of the boundary condition (3.1) was tested by modifying the
coefficient A in (3.1). A complete study has not yet been made, but
decreasing A by a factor of 2 (outflow boundary at 2.067 feet) increased the
number of time steps required for convergence by 560 and made only small
changes in the final steady state.

The quality of the steady state was assessed both by comparing solutions
generated by the two boundary conditions and by comparing the solutions

generated by the two positions of the outflow boundary. We first consider the


-24-

Table I. Non-dimensionalized U and V obtained from (1.1) and (3.1)

(y = .00088 ft.)

x U/U,, (1.1) U/U,, (3.1) v/v, (1.1) (W/u),, (3-1)
grid point
ll «1234 1234 «331 x 1074 331 x 1074
15 21144 +1144 +263 x 1074 .261 x 1074
19 -1071 .1072 218 x 1074 213 x 1074
23 .1009 +1014 «194 x 1074 174 x 1074
27 09526 09680 «201 x 1074 131 x 1074
31 .08978 09429 .231 x 1074 0821 x 1074

Table II, Non-dimensionalized U obtained from (1.1) and (3.1) with
outflow boundary at x = 2.07 ft. compared with values obtained
with outflow boundary 3 ft. from leading edge (y = .00088 ft).

x u/U, (1.1) u/U, (3.1) U/U,, outflow at 3 ft.
grid point from leading edge
5 1420 21421 1421
7 «1348 +1350 1349
9 -1286 . 1289 -1287
ll «1232 1238 1234
13 1184 21195 «1186
15 1139 1161 21144

17 +1106 «1143 +1106





-1ll-

£1 g ¢2)glV cam a - gf) - a2 \p =
( 09 YT ge )ss + ( +i y-1 ce) ong ~ ofl 7-1 “o}8 0.
(2.35)
The boundary conditions for (2.35) are, to zeroth order in 8:
= gt =
8 (0) = 896) 0 (2.36)
1 1
g(t) = 0% ay go(1) (2.37)
\, J
eee = Dy *h "
8 (1) io y-1 ory go(). (2.38)

In equations (2.35) - (2.38), 0 = Gps. For high Reynolds number flow,
say Rs > 30, € is a small number, and we can ask ourselves how does the
eigenvalue % behave as a function of ce ase diminishes and approaches
zero.

There are three possibilities:

0
lin = = yo .
£70 o constant

Another way to pose this is to say that wy ~ e* as ¢ + 0 and we ask whether
a is smaller, greater than, or equal to zero. If a < 0, then we have in
(2.35) a singular perturbation problem and one can show that, considering the
boundary conditions (2.36) - (2.38), the only way to match the inner and outer

expansions is to have the trivial solution. If a > 0, we have a regular


-13-

the graph) the |1/A9| curve deviates from the one shown by less than 5%.
Thus, we conclude that for moderate and high Reynolds numbers the very long
wave (8 * 0) decay rate is basically a universal function of o = (y-1)M2.
The correction, Buys to the decay is found by considering the first-order
problem, (2.32) - (2.34), with its attendant boundary conditions. Again, one
can extract from this system a single fourth-order differential equation for

& = of), After some manipulations the 8)~problem may be stated as follows:

Y 2\ IV 2y 2 Y 2

-—— —_ —)|e" - — -— — =

(% 7-1 oe )a + ( +a y-1 o " }e ongy (' 1 =} R (2.40)
0.

where

= we(o, R +R), (2.41)
R, = #2 - Qieg" - 20 (2.42)
1 2 2 8. 80 0 80 .
YM, w
0
(yt)
_ ite eecye _ ACytl) oon 2te 02 _ 20
Ro az (089 7 NBG = 8G yo Ng: (2443)
YM, 4%

The boundary conditions, to first-order in 8, are found to be
g,(0) = g}(0) = 0 (2.44)
Y, Y, Y, +1
- - 2(L) 7/2 = 2 eb
BCL) - do "(SF)" og 8) = 4G) 2 (u, + eg) (2.45)

ai) ~ to Ye (te ) Dw, ey) = 1) 2 (a, +P yep, (2.46)


NON-REFLECTING BOUNDARY CONDITIONS

FOR THE COMPRESSIBLE NAVIER-STOKES EQUATIONS

Saul Abarbanel
Tel-Aviv University

Alvin Bayliss
Exxon Corporate Research Science Laboratories

Liviu Lustman
Institute for Computer Applications in Science and Engineering

Abstract
A small perturbation analysis, in the long wavelength regime, is used to
obtain the downstream boundary condition for the pressure for the flow over a
flat plate. The methodology is extendable to other geometries. Numerical
results for high Reynolds number laminar flows show great improvement in

convergence rate to steady state as well as in the quality of the results.

This work was supported in part by the Air Force Office of Scientific
Research (NAN), AFSC, United States Air Force under Grant AFOSR-82-0136E while
the first author was at M.I.T.; and in part by NASA Contracts No. NAS1-17070
and NASI1-18107 while the authors were in residence at ICASE, NASA Langley

Research Center, Hampton, VA 23665-5225.
N8O-90 893%



(1)

[2]

[3]

[4]

(5]

[6]

[7]

-22-

REFERENCES

MacCormack, R. W., in Lecture Notes in Physics, Vol. 8, p. 151, Springer-

Verlag, New York/Berlin, 1971.

Olson, L. E., P. Re McGowan, and R. W. MacCormack, Numerical Solution of
the Time-dependent Compressible Navier-Stokes Equations in Inlet Regions,

NASA TM-X-62338, (March 1974).

Shang, J. S. and W. L. Hankey, Jr., Numerical Solution of the Navier-
Stokes Equations for Supersonic Turbulent Flow over a Compression Ramp,

AIAA Paper 75-4 (1975).

MacCormack, R. W. and B. S. Baldwin, A Numerical Method for Solving the
Navier-Stokes Equations with Application to Shock-Boundary Layer

Interaction, AIAA Paper 75-1 (1975).

Rudy, D. and J. C. Strikwerda, A Non-Reflecting Boundary Condition for

Subsonic Navier-Stokes Calculations, J. Comp. Phys., 36, (1980), 55-70.

Rudy, D. and J. C. Strikwerda, Boundary Conditions for Subsonic

Compressible Navier-Stokes Calculations, Computers and Fluids, 9 (1981),
327-338.

Bayliss, A., P. Parikh, L. Maestrello, and E. Turkel, A Fourth-Order
Scheme for the Unsteady Compressible Navier-Stokes Equations, AIAA Paper

85-1694 (1985).


3 1176 01358 3225

ean

‘


-25-

Table III. Non-dimensionalized V obtained from (1.1) and (3.1) with
outflow boundary at x = 2.07 ft. compared with values obtained
with outflow boundary 3 ft. from leading edge (y = .00088 ft).

x V/v,, (1.1) v/v, (3-1) v/v, outflow at 3 ft.
grid point from leading edge

5 1542 x 10-4 1538 x 1074 1538 x 1074
7 1448 x 1074 w41 x 10-4 1443 x 104
9 .388 x 1074 6374 x 1074 .380 x 1074
li .347 x 1074 4.320 x 1074 .331 x 1074
13 1325 x 1074 1268 x 1074 .293 x 1074
15 312 x 1074 .202 x 1074 1263 x 1074
17 2302 x 1074 .165 x 1074 .238 x 107*

1

Sval

9

2.6

2.0

3 (0 <e< .0053)

10° < R, <@
J L (y=1.4)

1.0



where r= r(M; w,8) is given by (2.19).

The system (2.13) - (2.15) with the boundary conditions (2.20) - (2.24)
constitutes the eigenvalue problem from which, in principle, we could extract
the dispersion relation w = w(f8) for given Mach and Reynolds numbers
(M, and ely,

As mentioned in the introductory section, we are interested in very long
waves, i.e., the small 8. We, therefore, propose expanding the eigenfunc~
tions G,; (i = 1,2,3) in the form
(0)

G, = G

(1) 2 (2)
7 + BGS" +B GE

toes (i = 1,2,3) (2.25)
where
of = [6,(n30,8,€)] 09° (2.26)
We also expand the dimensionless "frequency" w (which for small 8 will
turn out to be the decay rate) in a similar manner:
w= oy + Bu, + or (2.27)

where

w= [wl] gos (2.28)

Using the expansions (2.25) and (2.27) in (2.13) - (2.15), we get the

following zeroth and first-order problems in 8:


-12-

perturbation problem, and again the only possible solution is the trivial one
of Btn) +0 as e+ 0. We thus remain with the third possibility of
Wg ~e for small ec. Let this proportionality constant be e/u9 = oe

Then (2.35) is rewritten as:

2

o

yo 2\ Iv 2y e2 ¥
-_ i—_—— —yJ|er - ad -—_— =
(a y-1 © )ss + ( +i y-1 ° a ongy o(1 ya )éo 0. (2.39)

a

oN

It is assumed, and will be verified, that for reasonable Mach numbers, say
below hypersonic speeds, Xo = 0(1).

We can now state the "8 = 0 problem" more concisely: given o, ©, and
Y (i.e., the Mach number, Reynolds number, and the gas constant) solve the
differential equation (2.39), subject to boundary conditions (2.36) - (2.38),
for the eigenvalues ro and the eigenfunction & y(n) = of, This
eigenvalue problem was solved numerically for a range of values 0 < « < .05,
0<o< 20 (for laminar flow this corresponds roughly to Reynolds
number, 50< R, <* and Mach number, 0 <M, <7 for y = 1.4). For this
range of the parameters it is found that Xo is purely imaginary negative--
this means that % is purely imaginary positive. Thus, Q vepresents the
decay rate of very long waves--to be corrected for finite but long waves by
the term Bw,. Figure 1 shows a plot of [1/9 | versus o; for
O0<e< .0053 «03 <R, <<), Xo is completely insensitive to changes in the
Reynolds number and indeed may be taken to depend on the Mach number only.
As o +0, [1A] takes on the value r/4 which can be derived
analytically.

It is seen from Figure 1 that, for this range of parameters, [a is

ol
indeed of order unity, .4 <% Ido w5. For R, as low as 40 (not shown on