NYO-1480-167

Courant

r-

.,.,.

WEW YORK UNIVERSITY

c

WURAMT

Institute of

INSTITUTE- LIBRARY.

^^

Mathematical Sciences
AEG Computing and Applied Mathematics Center

Numerical Solution of Nonlinear

Boundary Value Problems Using Reflection

A. D. Snider

AEC Research and Development Report

Mathematics and Computing April 1971

New

York University

--

co^--TS.

AEC Computing and Applied Mathematics Center Courant Institute of Mathematical Sciences New York University

Mathematics and Computing

NY0-l480-l67

NUMERICAL SOLUTION OF NONLINEAR BOUNDARY VALUE

PROBLEMS USING REFLECTION

A.

D.

Snider

Contract No. AT(30-1 )-l480

UNCLASSIFIED

NEW YORK UNIVERSITY COURANT INSTITUTE- LIBRARY

TABLE OF CONTENTS

Part

Introduction

1
6
6

Part

2.

Reflection

A. Formulation of Free Surface Problems

B. Parametrization by a Conformal Transformation
C. D.

10

Calculation of Conformal Maps

Reflection
Implications of the Reflection Rules
and Convergence

13 16

E.

21
25

F. Existence, Uniqueness, G.

Summary
3.

27
30 30

Part

The Vena Contracta

A. The Physical Problem and Parameters

B. The Conformal Map
C. Application of the Reflection Laws to the X-axis

34

and the Wall

D.

35 36

The Reflection Law at the Free Surface

E. The Separation Point F.

49 54

The Line of Truncation

G. The Singularities at Infinity

H. The Interior Points
I.

57 64
66

The Fortran Program

4.

Part

Results

68
69

A. Dependence on Convergence Factors B. Accuracy of the Solutions

C.

71
73

Rate of Convergence of the Iterations

D. Estimates of the

Contraction Coefficients
-111-

75

Bibliography
Figures

77

78
I.

Appendix

Listing of the Fortran Program

86

Appendix II.

Dimensional Perturbation

98

-IV-

ABSTRACT:

Recently the solution of some nonlinear free

boundary problems has been effected niimerically by

incorporating a
change of independent variables through

a conformal transformation, thus simplifying the domain of

the solution (cf. E. Bloch,

[1]).

One then solves the

transformed differential equation in the simpler domain,

simultaneously determining the transformation itself, and

the result is regarded as a parametrized form of the desired

solution.

There is some question as to what is the best

way to handle the boundary conditions in such an approach

The aforementioned report employed a steepest descent

procedure to produce difference equations at the boundary

which uniquely determined the solution, but led to a rather

slow iterative scheme.
The present paper discusses the results of using
the reflection property of solutions of elliptic partial

differential equations [cf.

8]

to determine the boundary

conditions to the transformed difference equations.

The

basic idea is analyzed theoretically, and we demonstrate

its applicability to a special class of two-dimensional

problems.

The procedure is then applied to a simple

-V-

plasma containment problem and to the two-dimensional and

the three-dimensional axially symmetric vena contracta

models.

The results are compared with those of Bloch,

and it is seen that the reflection scheme requires about

one-tenth as many iterations for convergence as the method

of steepest descent.

New computations for the

contraction

coefficient are presented.

-VI-

Part

1.

INTRODUCTION
The advent of the modern computer has stimulated

the finite difference technique of finding approximations

to solutions of partial differential equations.

Nonetheless

the enormous number of calculations required by this

technique strains the capabilities of even the biggest

machines, necessitating studies seeking more efficient ways
of formulating and solving the equations
.

In the case of

elliptic boundary value problems one is usually confronted

with the task of solving a large system of algebraic equations expressing the finite difference analogues of the

partial differential equation and the boundary conditions.

The solution of this system is achieved by some iterative

procedure.

Thus the effectiveness and feasibility of

the

a finite difference scheme are indicated, by first,

rate at which the iterates of the trial solutions of the

difference equations converge to the actual solution,

and second, the degree to which this solution of the

difference equations approximates the solution of the

differential equation.

The latter is a measure of the

accuracy of the finite difference approximation.

The accuracy and rate of convergence of many

linear systems can be estimated theoretically [2,4,6,10].

The appearance of nonlinearities
,

however, makes these

In fact, in the case

estimates prohibitively difficult.

-1-

of nonlinear boundary conditions

particularly those

involving free boundaries, it may not be clear how to

formulate the finite difference approximation.
In this

paper we describe a method for handling some of these

nonlinear examples

Applying the technique to a free

bounaary problem, we are also able to establish the

accuracy and convergence of our scheme by experimenting with tne number of mesh points.
In the past,

free boundary problems have been

attacked by iterative procedures which guess the location

of the boundary, solve the rest of the equations, and then

improve the guess

[10]

A more accurate and faster scheme

He incorporates a change

was recently reported by Bloch [1].

of independent variables via a conformal transformation,

simplifying the domain of the solution.

Then he solves the

transformea differential equation in the new (known) domain,

simultaneously solving for the transformation itself.

result is regarded as a parametrized
solution.

The

form of the desired

The boundary conditions for the transformation

are derived from a steepest descent argioment, supplemented

by equations serving to establish the free boundary

constraint, which in his case expresses continuity of

pressure across a fluid-air interface.

In this paper we propose an improvement of this

procedure whereby the boundary conditions of the

transforma-

tion are obtained by analytic continuation, resulting in

-2-

reflection rules (cf.

[8])

replacing the steepest descent

The

and constant pressure equations of the above method.

scheme is just as accurate as Bloch's, but the iterations

converge almost ten times faster.

Furthermore, it is

directly applicable to any problem in which the solution

can be continued across the boundary in some prescribed

manner as described by H. Lewy in

[8]

for example.

This

report purports to demonstrate that the reflection technique

is superior to the method of steepest descent in handling

these problems
In Part
2

of this paper we describe the theory behind

It is introduced in the context of a

the reflection scheme.

general free surface problem of fluid mechanics, but the

actual procedure for obtaining the conformal map is presented

without such motivation.

Proofs of certain results are

given when available, and evidence of other expected

consequences is offered so that a heuristic understanding

of the scheme's implications and its limitations is obtained. In fact, the whole theory can be neatly capsulized, and

this is done in a summary at the end of the section. The application of the procedure in solving the

vena contracta problem is given in Part

3.

The physical

situation is the following:

one has an incompressible

inviscid fluid under high pressure confined behind a plane

wall, and the wall has an aperture through which the fluid

escapes as a jet.

In the two-dimensional model the

-3-

aperture is an infinitely long slit; in the three-dimensional

proDlem with axial symmetry, the aperture is a circular hole.

The boundary of the jet is of course unknown beforehand,
so we have to deal with a free boundary problem.

To determine the vena contracta one must modify

tne general procedure of Part
2

somewhat, so additional study

of tne theoretical aspects is necessary.

Furthermore, this

problem is full of complications not related to reflection,

and each of these must be dealt with.
In particular, the

behavior of the flow near the edge of the aperture is quite

singular, and we develop some special techniques for handling
this point.

Knowing the location of this edge, one can

evaluate the contraction coefficient, which is the ratio

of the area of the jet at infinity to the area of the

aperture.

Tnus we are able to present a new, accurate method

The details of all these

for computing this number.

considerations are presented in Part 3, which concludes

with a description of a Fortran program (listed in Appendix II)
for the solution
of the problem.

Part

describes the results of running this program

The accuracy

on New York University's CDC 6600 computer.

and efficiency are reported, and it is shown that the method converges about ten times faster than the steepest descent

procedure mentioned above.

The new calculations of the

contraction coefficient produced by the reflection program

are also presented in this section. They agree extremely well

with those of Bloch

[1]

-4-

The evidence offered herein clearly demonstrates

the superiority of this scheme in handling the vena contracta

model, and we are led to conclude that reflection is the

appropriate method for solving a wider class of free

boundary problems which arise in fluid dynamics and plasma

physics.

Furthermore, the results tend to encourage

the use of reflection rules whenever they are applicable in solving more general nonlinear elliptic boundary value

problems in two independent variables

As a final note, we present in Appendix
I

re-evaluation of a calculation of the axisymmetric contraction coefficient by a perturbational technique

[3]
,

which

motivated the present study.

-5-

Part
A.

2.

REFLECTION

Formulation of Free Surface Problems

One of tne principal sources of nonlinear boundary

value problems for which the reflection scheme is suitable

is

the free surface problem in fluid dynamics

so a survey

of the known results in this area will serve as our

introduction to the

method.

One must solve a partial differential equation for

the stream function
ip

in a region, called the flow region,

occupied by the fluid.

This region is bounded partly by

fixed walls and obstacles and partly by a constant-pressure

medium, such as the atmosphere.
The locations of the walls

and obstacles are known data, but the location of the free

surface must be determined as part of the solution.

value of
ip

The

is given on all boundaries,

and some other

condition (usually expressing the fact that the pressure

in the fluid must match the atmospheric pressure)
is given

on tne free surface.

Thus one is presented with a differential

equation to solve in an unknown region with, usually, nonlinear

boundary conditions.
Specifically, in two-dimensional potential flow
problems, the components of velocity in the (x,y) plane
are obtained from
i>

through the equations

,

= rr^ dy

= - -K^ dX

-6-

The equation

for

ij;

then expresses the condition that the

flow is irrotational:
-

curl V =

A4;

The flow on the fixed and free boundaries is parallel

to the boundaries,

so
Tp

= constant

there.

On a free boundary, the pressure p is constant.

Since Bernoulli's law relates pressure and velocity [2],

p +

p-

= constant

the velocity is also constant there.

Thus on the free

boundaries where the velocity is given by the normal

derivative
3i|i/3n,

we must have
=

'
r.

3n

constant

Of course the values of all these constants must be chosen

to satisfy conditions at the boundaries, at infinity, etc. In three-dimensional problems with axial symmetry,

we let x denote distance along the axis of symmetry and

y denote radial distance from this axis.

The velocity is

derived from the stream function

V
=

i)

according to the equations

i-|l
y dy

= -

iM y dx

If the flow is irrotational.

-7-

curl

v=--(Ai|j--|^) =0 y
y

and therefore
A4;

y 9y

^ |i =

Again, the flow on the free and fixed boundaries is

tangential, so on these Doundaries

= constant

ip

Bernoulli's law dictates constancy of speed on the free

boundaries as before.
1

This is now expressed by

d\b
T-"-

y dn

= constant

It has been shown [5,7]

that these problems have

unique solutions
are analytic.

and furthermore that the free boundaries

This latter fact will be helpful in

implementing the reflection scheme.

An example of this kind of problem, which will be treated in our paper, is the vena contracta.

Water is

kept under pressure in an infinite tank bounded by the

y-axis
wall.
,

and allowed to escape through an opening in the

The opening is a slit in the two-dimensional case,

and a hole in the three-dimensional axisymmetric case.

If we neglect gravity, the water escapes through the

aperture in a stream, which eventually looks like a tube

of fluid in uniform motion
(see Fig.
1)
.

The boundary

-8-

of the stream is, of course, unknown, but the pressure

at this boundary must equal atmospheric pressure. This

condition determines a unique solution.

As we mentioned
(x,y)

above, the boundary is an analytic curve in the

plane,

The two-dimensional problem can be solved by the

hodograph method.

The fiinction

^p

is harmonic and can be

c})

combined with its harmonic conjugate

fxinction
r,

to give an analytic

(i>

+ i^

The velocity is then given by

calculated it to be .59135 + .00004.

another evaluation below.

We will present

B.

Parametrization by a Conformal Transformation

The solution of free boundary problems can often

be aided by a parametrization via a conformal map.

is trying to solve an equation for
\p

One

in terms of x and y

in an unknown region.

Now if we regard the region as a

z =

domain in the complex plane of

x + iy,

and if it is

simply connected and not the whole plane, then the Riemann

mapping theorem tells us that this domain can be regarded

as the image,

under an analytic function

= z (w)

of

something specific, like a rectangle, in the w-plane.

Furthermore we can specify the boundary correspondences

of three points.
If we have this map, we can parametrize
ij;(x,y)

the equation for

to yield an equation for

i|;

as a

function of u and v, where w = u+ iv; this equation now

has to be solved in a simple, known domain.

The parametrization of the

ijj

equation in the case

discussed above is aided by certain identities involving

the gradient operator V and the Laplacian
A;

these

identities are immediate consequences of the Cauchy-Riemann

equations for the analytic fxinction

z (w)

If V

denotes

the gradient operator with respect to x and y while V

denotes this operator with respect to u and v, then for

-10-

any real twice-dif ferentiable functions f(x,y)

and g(x,y),

we have V^
f (x(u,v) ,y (u,v)
)

V^g(x(u,v) ,y (u,v)

and
A^

f(x(u,v),y(u,v))

1^1

A^ f(x,y)

Therefore, for two-dimensional potential flow

problems, the equation
A^
'\>U,y)

becomes
A

<J;(u,v) ' ^

The boundary values of ^ just become transferred to the

images of those boundaries.

To facilitate studying the

constant-speed condition, let us introduce the symbols

n
z'
,

and t

denote normal and tangential ^ w to

z

directions in the
_
z

and w planes.

Then by the chain rule,

9i|j

ij;

w w
9n
z

j)

w
w
9n
z

8n

9n

9t

But

9ip/9t
i|j

on the pre-image of the free boundary So the constant speed

because

is constant there.

condition becomes
9iij

V-'-

9n

\1;

-"r-

dn

TT

9n

w = constant.
.
.

-11-

since

z (w)

is conformal and maps boundaries to boundaries,

the normal to the bounaary is preserved and we can write

3n
3n
'^

9n
z

3n

'

3n

'

provided the positive directions along the normals are

chosen appropriately; thus we can rewrite our condition as

^r^ / Art 3n w
'

i '

l-:r^ ' Ar\

9n

constant

on the free boundary.

The axially symmetric equation for

i|j

is

parametrized

just as easily when it is written in the form

V

y
5

'

4^

AW;
z^

_ =
y

It transforms to
V y

i^

w^

]h

y
ip

in the w-plane
as before,

The boundary values of

are carried over

ana the constant speed condition becomes

y 3n
-^

- /
'

TT

3n

'

= constant.

w
ip

So now the problem for

is simpler,

given the

analytic map; one has an equation with boundary conditions

similar in form to the original, but these are to be solved in a known region.
The fly in the ointment here is that one must also
find the conformal map.

And when there is an unknown boundary.

-12-

as in free surface problems,

the mapping problem is

ip ,

coupled with the problem for

through the free boundary

condition (as in \^\ / \^\ = 1); so one cannot, in

general, solve the equation and find the mapping separately.

Therefore we shall examine the conformal mapping problem

in detail.

C.

Calculation of Conformal Maps

Let
z (w)

be analytic; then, of course, x and y

are harmonic functions of u and v:

Ax =
,

Ay =

Now reversing our point of view, if we wish to find a

conformal map, we can start by regarding these as

differential equations to be solved for x and y;

but then the boundary conditions must be determined to

fix the solution.

First of all, the image

= x +i y of a point on

the boundary must lie on the curve bounding the image

region.

This is a rather subtle requirement, however,

because it does not specify x and y individually; at most

it

gives one in terms of the other, and in the case of

an unknown free boundary, it tells us nothing. More

conditions are needed.

In
[1],

Bloch utilized the theory of Plateau's

As

problem to get the necessary boundary conditions.

-13-

described in

[2]

if one has continuously dif ferentiable

functions x(u,v) and y(u,v) defined on a rectangle in the

(u,v)

plane and mapping this rectangle onto a domain bounded

by a closed curve, and if one minimizes the Dirichlet-

Douglas integral
((Vx)^ +
I

(Vy)^)

du dv

R
subject to a three-point normalization (i.e., the images
of three points on the boundary of the rectangle are fixed)

then the resulting function x+iy of the complex variable

u+iv is unique and analytic.

To minimize the integral

under any auxiliary conditions, Bloch shows, one must

have Ax = Ay =

and also on the boundary of the rectangle

(x,y)

the normal derivative of the vector

must be
One can express

perpendicular to its tangential derivative.

this condition as

x^x

t n

+ y.y
-^

t-'

where n and

denote normal and tangential derivatives

The left-hand side is the dot product of
,y
)

respectively.
the vectors
(x

and (x ,y^)

In the case where the

rectangle in the u,v-plane has sides parallel to the axes,

the orthogonality condition asserts that

XX u V
on all sides.

+ y y u-' V
-^

=0

-14-

To understand exactly what this last condition

says about the solutions to the equations Ax = Ay =

one must investigate it further.

Following Bloch, we
is harmonic if x and y

note that the function x x +y y

are, so if it vanishes on the boundary then it vanishes

inside the rectangle as well.

2 1,2^2 -X u ^(x 2 V +y 'v

Its harmonic conjugate is

2.
)

- Y 'u

= -^
2

Ij,
I

as direct calculation shows

so this latter function must

be constant in the rectangle. Furthermore, if this

constant is zero, then the map x+ iy is conformal

because it satisfies the Cauchy-Riemann equations.
In summary, if the Dirichlet-Douglas integral is

minimized subject to a three-point constraint, the map is

analytic and the constant His zero; and if it is minimized
subject to additional constraints but with H =
,

the

map is still analytic.

Bloch actually solves the problem by enforcing
a four-point condition specifying the images of the corners
of the rectangle and iterating his answer by a steepest

descent correction according to the following:

1.

moving (x,y) along the boundary in such a way

as to drive x x +y y

u V

-^

u-*

to zero,

2.

moving the free boundary itself in such a way

as to drive the pressure to its correct value,

3.

moving the (prescribed) position of two of the

-15-

corners in such a way as to drive H to zero.

This scheme is successful, but it is our ooal here to

establish that reflection is equally successful and much

more efficient.

D.

Reflection
To describe the use of reflection in seeking a

conformal map numerically, let us consider a rectangle

in the w-plane. It is to be mapped onto a region in the
2)
.

z-plane bounded by real analytic curves (Fig.

To

find the map, we observe that x and y are harmonic functions

of u and v, and we use this in a finite difference scheme.

Let us first set up a mesh in the rectangle (Fig.

If we identify the points by integer pairs

3)

(i,j), enumerating

them in the u and v directions

respectively, then for

harmonic functions we have the approximation

X.

1,3

4 T

(x.,,

.+x.
D

.+x.
D

1+1,

1-1,

i,D+l

.,,+x.

i,D-l

,)

+ 0(h

2
)

where h is the mesh size, assumed uniform.

a similar equation holds for y.

Of course,

This is the difference

.

equation approximating the equation Ax =

This finite difference analogue of the underlying

differential equation can be written at every interior mesh

point of our rectangular grid.
In standard treatments of

boundary value problems, e.g. the Dirichlet problem, another

-16-

difference equation approximating the boundary condition

is written at the mesh points on the edge of the rectangle.

However, in nonlinear examples such as the conformal mapping

of the rectangle onto some arbitrary region this may not

be feasible.

To circumvent the difficulties we introduce

new mesh points just outside the rectangle and use them
to write this same difference analogue of the differential

equation at the boundary, as well as the interior, mesh

points.

Then the problem becomes finding reflection

rules, based on the boundary conditions, which will

determine the values of the solution at the reflected

mesh points.

This procedure is familiar in treating the

Neumann problem, and it usually leads to a convergent

iterative scheme for finding the solution even when

nonlinearities occur.
Clearly the technique is applicable, in theory
at least, to any problem in which the solution can be

continued across the boundary according to known rules.

In
[8]
,

Lewy demonstrates this property for a class of

analytic elliptic equations, using the appropriate Riemann

functions.
In our case, where a straight line in the
z

w-plane is mapped onto an analytic curve in the

-plane,

the property of analytic continuation is a well known

result of the theory of conformal mapping

[9]

So now we form additional mesh points outside the

rectangle which are reflections, through the sides, of

-17-

those mesh points adjacent to the sides, and x+iy can be

defined at these reflected points (at least, for small

enough
h)

(Fig.

2),

Thus we can write the numerical analogue

for Laplace's equation at each boundary point, even at the

corners; each one has a left and right and an upper and

lower neighbor.
This leaves us with the problem of writing

equations for the values of x and y at the reflected points,

that is, the formula for the reflection rule.
To derive this

rule, we must consider more carefully the exact description

of the analytic curves bounding the image region in the

z-plane
Let us say the real-analytic curve
T

is one of

the bounding curves in the z-plane, and its pre-image is the bottom of the rectangle, v =
for
r

(Fig.

2).

The equation

is of the form

g(x,y)

where g can be expanded locally in a real power series.

Making the substitutions

z+z z-z

(where a bar denotes complex conjugation) we get the

equation for

in the form

F(z,z)

where F(z, ,z) is complex-analytic in

z,

and z_ because

F has a power series expansion, namely, that derived from g.

-18-

Now if w_ is an interior point in the w plane and w

is the
0,

reflected image of w_, reflected through the line v =

if z_ = z(w_)

and

and

z (w

then we assert that

.

F(z^,z_) =

This is the basic formula which gives values of

at the

reflected points in terms of the values at interior points.

Before we prove this, let us examine two special
cases
Case
1.
r

is a straight line,

say x = 1.
0,

Rewriting this as x-1 =

then as

we get
z

+z
1

F(z^,Z2) =

-V^-

We know the rule for reflecting in this case; it

is
(z

-1)

= -(z_-l).

This agrees with F(z, ,z) =

2

Case

2.

is a circle,

say x

+ y

=1.
and
.

This is written as zz-1 =

F(Zj^,Z2)

= z^Z2 - 1

The rule for reflection is

1

z_

again agreeing with F(z ,z_) =

0.

-19-

Of course, if F(z,z)

is to specify a curve,

then Re F(z,z) and Im F(z,z) are not independent functions;

if they were,

the (complex)

condition F(z,z) =

would

define, as a locus, the intersection of two different

curve:;,

i.e., certain isolated points.

In the two cases

considered above, Im F(z,z) was identically zero, as it

will always be if F is derived from g in the manner
described.

Now we prove the rule for reflection.

z

Define

and w as before.
.

Theorem

If

z (w)

is analytic and if it is defined on both

sides of the line v = 0, mapping this line onto the analytic

curve F(z,z) = 0, then F(z ,z_)

- 0.

Proof

Noticing that w

and w_ are related by w_ = w

we define the function

f(w)

= F(z(w) ,zCw))
z (w)
.

where

z (w)

denotes the conjugate of

f (w)

When w = w

w = w_ and

is just F(z

z_)

We know F is analytic

in its first and second variables.

Since

z (w)

is also

analytic,

z (w)

is analytic and thus f is analytic in w,

by composition.

When w lies on the bottom of the rectangle,

z

w = u = w, then
f(w)

lies on T, so
= F(z,z)

= F(z{u),i(u))

Therefore

f (w)

an analytic function of w, vanishes on

-20-

the line v = 0; hence it is identically zero.

Thus

F(z^,z_) = f(w^) =

Q.E.D.

From the above proof, we can readily see that the

same rule holds for reflection through the other sides of the

rectangle also, i.e. F(z^,z_) =

is the general rule for

analytic continuation

Now we have a feasible scheme for calculating

conformal maps.
We write the numerical equivalent of

Laplace's equation at each interior and boundary point,

and we write F(z
,

z_)

for each reflected point.

The next section discusses the consequence

procedure

E.

Implications of the Reflection Rules

Now we must reverse our point of view in this sense;

we set up

Baplace's difference equation at each interior

and boiindary point of the rectangle, and we write the

reflection rule for each of the reflected points

can we

What

say about the solution to this system?

For the moment let us assume that the solution

exists and is unique for every sufficiently small mesh size h,

and that the solution converges to a function z(u,v)
= x+iy,

where x and y are harmonic, as h goes to zero.

These

assumptions will be discussed in the next section.

the properties of
s

What are

on the boundary?

-21-

Let us examine a particular point w on the boundary

of the rectangle, together with its inner neighbor w_ and
its reflected neighbor w the three points
z

(see Fig.
)

2)

As h goes to zero,
z_ =
z (w_)
,

= z (w

z-

= z(Wq),
z,
.

must

all converge to the same point

But F(z ,z_) = 0, so

F(z,,z^) =

and

z,

lies on F.

Another way cf seeing this

is to notice that z

and z_ always lie on opposite sides

F
.

of

r,

so

z z

must lie on

We conclude that the limiting

function

maps the boundary of the rectangle onto the

boundary of the region in the z-plane.

Now it is our goal to prove that the normal and

tangential derivatives of
X X +y y = 0.
z

are perpendicular, i.e.,

To this end we must study the geometry of

the curve

First

v^e

assume that the function F(z,,Z2) was

derived from a real function g(x,y) as described earlier,

so that F(z,z)
is always real.
z

If
,

is given as a function

of a real variable s,

= z(s)

then

dF
-a

must be real.

Therefore
dF = F,z 5^
ds
is real.

Is

+ F-z
2

Since z(s)

is arbitrary,

we conclude

Re F, = Re F-, Im F, = - Im F_; in other words

Now if z(t) is a parametrization of the curve T,

F(z(t) ,z(t)
)

= 0.

Thus

-22-

g=

= F^z^+ F^i^ = F^z^^ F^i^ = 2 ReCF^z^)

= Re F^x

+ Im F2y^

Now

(x

,y.)

is tangent to the curve r, so

(Re Y^,

Im F2)

is

normal to the curve and (-Im F2,Re F2)

is again tangent to r.

If we associate, with each complex number, a vector

v/hose X and y components are the real and imaginary parts

respectively, we can say that F2(^f^) is normal to the

curve F(?,t) = 0.

We can get more information about

by applying
If
=

the analytic form of the implicit function theorem.

F(z,,Z2)

is analytic in both variables with F(a,b)

and F,(a,b) ^ 0, then F(z,,Z2) =

defines z^ as an analytic

function of z in a neighborhood of (a,b)

(z,-a)
= Cq + c,(z2-b)

+ C2(z2-b)

...

Furthermore, c =
We have F(z Let us suppose that we just proved
r>

and
,

c,
)

= -

^2

(a,b)
z

=0

defining

in terms of z_.
?

and z_ both approach

T;

as h
= 0.

0;

was on the curve

so F(^,C)

The theoren then tells us that we can write

F2(?,C)
(z.-s)

= -

(z -I)

- 2 + + c^(z_-C)

(CO

Using F,(5,5) = Y^it.iO and adding ^-z_ to both sides,

we ultimately have

-23-

'

As the mesh size h of the rectangle becomes smaller, ?-z_ ^+~^_ should give approximations the quantities -,and -rto the derivative of
-^

in the direction norraal to the side

z

of the rectangle, denoted

Thus if we divide the above

equations by 2h and take the limit, we get

z

"

^2

i]^

T [F^z n +F~z n 2 2

Since the coefficient of F^ is real, the direction of the

vector
i.e.

n'

if it is non-zero, is the same as that of F-, ' 2

T.

normal to the curve

(A case

where

is zero

will be mentioned in the next section.)

Of course, the derivative
z

of

along the side

of the rectangle will be parallel to the tangent to T, so we will have

^n^t
for the function z(u,v)

-^

^n^t =

along the boundary of the rectangle.

Reasoning as before for these harmonic functions x and y,

we know that the latter equation holds throughout the

rectangle and that x ^
2

2 2 2 -y +y -x u =His 'v -^u

constant there.

From what
that
z

v/e

have said so far we cannot conclude

is analytic, i.e.,

that H = 0.

For example, the

mapping x = 2u, y = v satisfies all the above requirements

-24-

and is not analytic.

To get an analytic map we must impose

some additional condition to drive H to zero; hopefully

the particular application will suggest such a condition.

F.

Existence, Uniqueness, and Convergence

We shall now discuss the existence and uniqueness

of the solution, and its convergence as h
-*
,

on the

basis of nmmerical experiments.

The method has been programmed for several different

cases of image regions bounded by four analytic curves.

The vena contracta described below leads to just such a

region.

The reflection rules were applied so as to map

each side of the rectangle onto one of the bounding curves.

This amounts to a four point normalization because each

corner must be carried into the intersections of these

curves.
The equations were iterated using overrelaxation

to accelerate convergence and in each case the iterates did

converge.

This indicates (but does not prove)

that the

solution of the four point normalization problem exists

and is unique, for if there were two solutions, one would

expect the iterates to oscillate between them.

By shifting

one of the curves one can get different values for H;

only one position makes H = 0, and then the finite difference

analogues of the Cauchy-Riemann equations are satisfied.

No proof of convergence of the method as h ^
has been developed, but the problem has similarities to the

-25-

Dirichlet and

Neijtmann

problems, whose finite difference

In the latter cases,

analogues have been analyzed [2,6,11].

the solutions of the difference equations approach those

of the differential equaton with, an error of the order

so we expect the same of our scheme also.

2 0(h),

Nimerical
)

experiments confirm that the convergence is 0(h

(see below)

Let us briefly examine the three point normalization.

If we choose to map a rectangle onto a domain bounded by

three analytic curves, we would naturally map two adjacent sides onto one of the curves, and we would not be specifying
the location of the included corner.
In Figure
4

we have

labeled the corner w

while its neighboring mesh points on

and w^, and on the vertical edge,

the top edge are denoted w,

w^ and w..

The neighboring interior point is called w_, and

its reflected images through the top and side are called

'

and w|

respectively.

Observe that the two reflections

z
,

of w_ are mapped into one and the same point,

by the

rule F(z

z(w_)

= 0, because the same reflection rule is If we denote the

applied to both sides of the rectangle.

image of each subscripted w by a
z

with the same subscript,

we write Laplace's difference equations

^1= i

(z_+Zc+2++23)

so 4z,-z^ = 4z2-z..

Therefore we can write

-26-

3z -4z,+z^ 3z -4z+z. 2 4 1 3 _ c c

2H

2H

The left and right hand sides are finite difference

approximations of order 0(h) to ^ and

8 z

at the corner,

so if we can assume that these converge to the derivatives,

we have

= z

there.
is
z

Thus x

= x

and y = y , so -'u -^v'

2 2 2 2 H = X +y -X -y V 'v u 'u

zero at the corner, hence zero evervwhere,

and therefore

is analytic.

We conclude that if the

reflection scheme converges for the three point normalized

mapping on a rectangle, the limiting function is analytic.

Notice that nowhere in the above proofs do we assume
that the image of boundary mesh points lie on the boundary of
the image region.

Only in the limit as h

-^

can we be sure

of this correspondence of boiuxdaries

G.

Summary
The results of the preceding section can be stated

concisely as follows.

Consider analytic boundary curves described by

equations of the form
F(z,z) =
.

We write Laplace's difference equation for x and y in a

region in the w-plane; the difference equations are vritten

for each interior and boundary point.

Then we write

F(z^,z_) =

-27-

for

tlie

reflected points.

This procedure constitutes the

reflection scheme for conformal mapping problems.

When the w-region is a rectangle and the image is

bounded by four analytic curves, each corresponding to a

side of the rectangle, our numerical experiments lead us
to conjecture the following theorem.

Theorem

The solutions

x,

and

y,

of the equations of the

reflection scheme converge to harmonic functions x(u,v)

and y(u,v)
as the mesh size h is diminished.
2
)

The rate

of convergence is 0(h

Furthermore, whenever the solutions and their

difference quotients converge, we have proved the following.

Theorem.

Under the convergence assumptions stated, the

z (w)

function

maps the sides of the rectangle onto the

z

boundary curves, and on these sides the derivatives of

in the normal and tangential directions are perpendicular

to each other.

As we mentioned, it follows immediately that

XX u v

+ y y u-' V
-^

=0

in the rectangle, and therefore X 2 + y 2 - x 2 - y 2 = constant V u u V

-'

Experimentation shows that the value of this

constant depends on the locations of the curves in the
z^plane, and the constant can be made zero by choosing

2 8-

the locations appropriately.

This condition will then

make

z (w)

analytic.

The above procedure seems to be most suitable for

applications of the method of reflection when the curves

in the z-plane are analytic and known.

However, as we

proceed in Part

it will become clear that neither of

these conditions is essential, and one can often adopt

the technique to less restrictive cases.

-29-

Part

3.

THE VENA CONTRACTA The reflection technique will now be employed to

solve the two- and three-dimensional vena contracta problems.

This section describes the details of the procedure.
A.

The Physical Problem and Parameters

As described earlier, for the physical situation in

the vena contracta, one has an infinite volume of incompressible

fluid under pressure behind an infinite plane wall; the wall

has an aperture which is either a slit (two dimensions) or
a

circular hole (three dimensions)

A jet of fluid escapes

from the aperture and its cross-section contracts as it moves

away from the wall; asymptotically, it becomes a tube of fluid
in uniform motion. In eitlier case, one needs only to examine the flow

in an (x,y) -plane.

For the infinite slit, one has transla-

tional symmetry in the z-direction; for the circular hole,

one has rotational symmetry about the x-axis

The motion is described by the stream function

'Jj(x,y),

which gives the x and y components of velocity by

V
V X
=

il / y^ 3y
'^

-^

= - Mi / 3x ^

y"^
-^

where m =

in two dimensions and m = 1 in three dimensions.

41

The equation determining

A^ -

is =

E |i
y 9y

-30-

Let

tiie

y-axis be the wall, with the jet issuing

to the right and the aperture extending from y = Yq to

y = -Y
.

Then the motion is symmetric about the x-axis

and we can merely consider the problem in the upper half

plane
Since the x-axis is a streamline, we can set
there arbitrarily.
one streamline, so
^p
i|^

The wall and the free surface also form

will be a constant,

iJj-,

there- The

fluid behind the v/all is under a pressure p^ at infinity,

and atmospheric pressure in front of the wall is a constant,
p.

The free streamline approaches a horizontal line, y = Y^,

asymptotically.
be a constant, U^

The speed of the flow on the surface must

,

by Bernoulli's law, so if n denotes the

direction normal to the streamline, we have

ii=
9n

in two dimensions, and

1
3'lJ

_ U

in three dimensions on the free surface.

These constants are related, of course.

The speed

U depends, via Bernoulli's law, on the difference p^- p

The value

i|^^

determines the rate at which fluid escapes in

the jet, so it depends on the other parameters also.

To obtain a consistent set of values, we observe the following,

-31-

First of all, the behavior of

jet is simple.
In (m+2)

^i

at infinity in the

dimensions.

for constants a and 3.

This follows from the differential

ip

equation and the fact that

Also, for large x

will ultimately not depend on x.

7n -

Thus the boundary conditions imply

a = Uq

= Uq/2

Since ^ approaches

i|^^

when

y = Y^

and x increases, we have

Uq^co

m =
ni

^
Therefore
ij^_

- 1

can be calculated from Y^ and (p^-p

Secondly, it is clear on physical grounds that

specifying Y_

p^, and p. completely determines the solution;

these are precisely the variables which the experimenter can

control.

Furthermore, only the difference p^-p^^ can be

relevant; an equal increase in pressure on both sides of

the wall cannot change the flow of an incompressible fluid.

-32-

Third, the geometry of the flow region, in particular

the location of the free surface, depends only on Yq.

Given

one flow with its associated Y^ and P<-Pq (determining Uq and

ip-)
,

we can find another flow for a different value of the

Poo~Pn
,

pressure drop
u'

^^"^
i|Jq

with

tlie

same Yq

by calculating

from (p1-Pq)

then
i|j

from Uq and the unchanged Y^, and

I

finally multiplying

by the constant

^Q/i>Q-

This flow

issues from the same aperture Y

and satisfies the pressure

condition appropriate to

^"^ because of the uniqueness p^-p^r 00

[J

theorem fV], it is the only one.

It utilizes the same flow

So, as we claimed,

region, hence the same free surface.

Yf.

alone determines the free surface, and, consequently, Y^.

Fourth, notice that a simple magnification of the

flow region for Y^ gives a flow region for a different Y(and hence the unique one)
,

so the non-dimensional ratio

Y^Yq
The

is independent of any of the parameters of the problem.

number iX^/Y^)^'^^

(m,

as before, is

or 1 in two and three

dimensions respectively)

is called the contraction coefficient,

and we shall determine this number in both the two- and three-

dimensional cases.
It follows that we can determine the flow region by

specifying Y^ instead of Y-; furthermore, this convention

will facilitate programming our solution, and we adopt the

following normalizations.
Y oo = 1 and U = 1
= It follows that ^n ^0
}
]

in both two and three dimensions.

"^~,

,-'
.

m= 1

p_ and p- are not needed

'^ "

for the solution (since we already have U_)

and Y^ must be

calculated as part of the problem.

-33-

The Conf ormal Map

As described in Part 1, the problen will be attacked

using a confoiT:.al map.

Theoretically, we would like to map an

The part of the

infinite rectangle onto the flow region.

rectangle extending to infinity would correspond to the jet,

the top of the rectangle would represent the free boundary,

and the bottom would represent the x-axis

The remaining

side would be mapped onto the fixed wall, and a simple pole
at the lower corner would map this corner to infinity in the

second quadrant. Figure

5,

These correspondences are depicted in

where ABE is the infinite rectangle.

In practice, we will approximate this situation

with a finite rectangle, truncating ABE with the vertical

side CD. The boundary correspondences will now be defined

precisely.
The corner A of the rectangle, at the origin of the

w-plane, will be mapped to infinity in the second quadrant.

This requires a pole in the transformation
z (w)
,

so near

w =

we have
z(w)

^ R/w
The

with a negative residue R which must be calculated.

side AB will be mapped onto the wall from y =
<

to y = Y.

The side BC is mapped onto the free boundary r.

So the

angle at the corner B is expanded from

Tr/2

to

it;

this is done

-34-

because near the corner

ip

has an expansion in (z z)

1/2 '
,

and the doubling of the angle which we will get keeps

the truncation error down to
(h );
2

see

[1]

for details.

The side AD is mapped onto the x-axis The image of the line CD is a curve

CD' which

becomes asymptotically a straight segment as the rectangle

is lengthened.

We truncate this behavior and ass\ame C^D*

is a straight line, but its location must still be determined

(Fig.
5)

The treatment of all the boundary conditions around

the rectangle is a complicated affair, and we will devote

the next few sections to this matter.

C.

Application of the Reflection Laws to the x-axis

and the Wall. The sides AD and AB are each to be mapped onto known

curves, so the reflection technique is appropriate here.

AD is mapped to the x axis, whose equation is z-z =

0,

Indicating images of interior points adjacent to AD by x_ and

y_, and their reflections by x

and y

(Fig.

5)

we derive

x^ = x_

and

y^ = -y_

AB is mapped into the y axis, and we have a similar

situation.

The reflection rule becomes

y^ = y_

and

x^ = -x_

-35-

Arguing as in Part

1,

we conclude that these equations

will demand that the boundaries do correspond and that

XX u V
along AD and AB

+ y y u-' V
-^

D.

The Reflection Law at the Free Surface BC is to be mapped onto the free boundary

r,

but the

situation is more complicated than that treated in Part 1.

The location of the free boundary is unknown, and we have

no simple equation like F(z,z) =

to work with.

However,

we can derive a r^^flection rule from the integral equation

for
ijj

in two dimensions, and we will modify it to make it

suitable in three dimensions.

Then we

v/ill

have to study

the implications of such a procedure all over again.

To obtain the integral equation for

i) ,

we again
2

employ the implicit function theorem as in Part

the equation F(z,z)
=
z

to rewrite

in the form = g(z)

where g is an analytic function.

With this formulation,

that the stream function

Garabedian has established in

[3]

for gener.ilized axially symmetric flow in (m+2)

is given by

dimensions

4;

(x,y;m)

= Im

2m
1

z
i

^:^_^^r"/2.._^,^NNm/2 (z-t)"'^^(z-g(t))

^1

^(_m^_m.
-36-

^.

(z-t)(z-g(t))
(z-t) (z-g(t))

)g.(^)l/2

^^^

where

is any point on the free surface T,

and

.j~

is

the hypergeometric function

~
^

r(i- ^)^ T

r(-

J)^

j=o r(j+i)^

The integrand on the right-hand side, considered as a function

of z, is essentially the Riemann function of the equation

8z 3z

2(z-z)

dz

for
^(x,y)
= 4^(^7^/ ^^f)
/

and this equation is just the complex forr. of

which is the usual equation for the stream function.

The integrand, considered as a function of t, is analytic
so the choice of
tlie

path of integration is immaterial.

^p

Finally, we add that the normalization for

in this
tlie

equation differs from ours, where

4/

= 1 on r, but

additive constant will drop out in our application.

The equation is particularly useful in two dimensions,

m =

it becomes
ij;(x,y)

= Im

g' (t)

dt

^1

In two dimensions,

i)

is harmonic. Thus it can be combined

-37-

with its hanT\onic conjugate

function
Z

(()(x,y)

to give an analytic

<t>

The above formula immediately

displays this property, since the right-hand side is

analytic, and
z
<t>

+ i ^ ^

g'(t)^/^ dt

^1

Differentiating, we see that

z

g(z)

(^)^ dz

^2

for some z-.

Staying with the two-dimensional case,

w_, w
,

we employ our usual notation:

and w

denote

interior, boundary, and reflected mesh points in the

rectangle, while z_, z-, and

are their images

(Fig.

5).

The free boundary equation says

Zq - g(zQ)

=
=

and the reflection law F(z

z_)

takes the form

.

z^ - g(z_)

=0

Siibstituting for g(z)

and subtracting, we have the reflection

rule for the two-dimensional free boundary problem.

z, ^+

- zZq = g(z g(zj)

- g(Zrt) g(zQ)

?'

dz

^0

It is now our task to approximate this rule so that it can

be used in our finite difference scheme.

-38-

We can integrate the right-hand side by employing

the power series development of
z,

in z, and that of

in w.
I

Expanding ?'(z) about

we have

and using the notation C'(Zq) = ^q,

By the chain rule.

^ d((j)+ij^l
1

^'

dw

dz/dw
Zq=z(Wq)

3v

dV

'

But since

({>

ip

on the top edge of the rectangle.

^0 = "^v/'o

Similarly we find
2

1
'0

^UV
^0

VV

'P

^v

z.

dz

"T77

where we use the Cauchy-Riemann equations to eliminate

Notice that since ^ is harmonic.
^

(j)

^vv

= -

^\,

^uu

and the latter is, again, zero along the top edge of the

rectangle.
to derive

Now we use the expansion for

in terms of w

-39-

"2
Z^ =

dzi
;;t

dw'w^

= -

^- ^0
i

r:

^0^ - ^ ., 2ih

^,,2. ^ + 0(h )

in terms of the mesh size h.

When all of this is substituted

It

into the integral, the appearance of z. cancels and we get

-h
z^ - Zq =
2
4i

+ih

3
iiitJj
-,

+ 0(h

From this approximate equation we formulate our

reflect:! on rule as

= z^ +

-h
2

ill

+ih ^ !z !ziHZ
\1)

The i/^-derivatives can be approximated by finite differences

accurately so that the truncation error on the right hand

side is of order h
2

higher

th.an

the left-hand side itself,

which is, of course, 0(h).

Heuristi cally this seems desirable

because we are using a difference equation which approaches

Laplace's equation to 0(h
2
)

and we should try to match this

accuracy in all our other formulae.

However, to renlly evaluate the effed of this rule,
we must studv its consequences as we did for the exact rules
in Part 2.

We assume we have a solution of Laplace's differ-

ence equation defined on the interior and boundary mesh points, At the points
z

which are reflected through

2
iJj

T,

we have

z_

-h
=

^v

+ih

3
li

^v uv

lii

^_~^o
(z
-z^^)
,

If we multiply by

divide by h

2
,

and take the real

part, we get

-40-

^^

-T
j-

= -

Notice that

r nand

are both difference approxiif we assume convergence

mations to

^ on

the boundary,

of the scheme as we did in Part 2, we can conclude that

the limiting form of the equation

_ Tz 9z ^
,2

or equivalently

</Nvl

- 1

is a consequence of the reflection rule.

There is another condition imposed by this rule,

corresponding to the imaginary part of the equation.

a little algebraic manipulation, we can derive
z

With

-2z+z_ z_-z
h
^ 72
r:

Im

ll;

^v^uv

it

Here we have a difference approximation to the second

~2

derivative

z/9v

so in the limit as h

->

the assumption

of convergence compels us to conclude

- Im z

vv

il;

^v ^uv

ill

as a second consequence of the reflection rule.

Before we go further, we should establish that

these consequences are consistent with the solution we
seek.

First, notice that if we have a conformal mapping,

ip

then the normal derivative of

becomes 84;/8n =

^ /\z.

and our constant speed condition says that

-41-

dn

This verifies the first consequence of the rule.

Second, the constant speed condition must hold all

along the top edge of the rectangle, so we can differentiate

with respect to u.

Expressing the condition as

z

V V

.2
lil'

we have
:

UV V
z

V UV

I'd)

^V UV

\b

If

z (w)

is analytic,

= -iz
2

and the left-hand side can

ultimate 7-Y be written as - Im

z

Im
ii

(z

Thus we have

^v^uv

ii

which is

just the second consequence of the rule.

We

might add that, had the reflection law not kept the h 3 term
with
^Ji'y.'v'

^^ would be imposing the condition

Im
z

=0

along

T,

which would be wrong.

Although we have shown that these two consequences

are consistent with our solution, they do not appear to be

as strong as

the ones generated by application of the exact

z

reflection laws; namely,

x^x^+ y^y
=
.

maps bovmdary to boundary and

In fact niomerical experiments on a very

simple plasma-containment problem indicated that the

application of the inexoct rules does not determine a

unique solution to the difference equations, as does the

-42-

exact formulation (with four point normalization)

Also some of the solutions calculated with the inexact
rule allowed the derivatives
u V
z

and

to become

non-orthoaonal ' , so that x x +y y was not zero. u-' V

-^

Further experimentation led to the suspicion that

the latter fault was the crucial defect of the procedure.

The difficulty is resolved by modifying the reflection rule

in such a way as to enforce the orthogonality of these

derivati\'es

One can accomplish this by adding an

.

appropriate function of the derivatives

The new term

should have the following properties:

(1)

it should cause
so as to make
z

to be shifted in a direction

and

more nearly perpendicular,

and the amount of the shift should be greater

when the non-orthogonality, as measured by

XX u V
(2)

+ y y , is greater, ^ u-' V '

-^

it should be zero when x x + y y is zero, so ' u V ^u-'v

that under this condition the term has no effect

and the two consequences of the reflection rule

will hold.
(3)

it must not cause instability, or

else the

iterates will still diverge.

The first two considerations immediately suggest that
the term be proportional to x x + y y
,

which, of course,
z

depends on the cosine of the angle between

and

-43-

To determine the direction in which

should be
z

shifted, we refer to Fig.

6,

where

and

denote

the images of three consecutive mesh points on the boundary,

z

and
z

denote interior and reflected points as usual,

and
z

Then

correspond to
z

ct and
z

^r-

respectively.
,

between When the angle ^

shown,
z

and

is less than tt/2

as

should be shifted in a direction 90 counter^

z

clockwise from the direction of

-z_, that is, in the

direction of i(z -z_).

than

When the enclosed angle is greater

in the opposite direction; the

correlated with that of x x + y y . sign change is exactly ^ uvu-^v Ultimately one sees that the simplest way of
introducing this shift into the reflection rule given above
is through a term of the form

^^

7t/2,

one shifts

X X +y y u V U-' V
-^

Here we have approximated the direction of

z-z_ because we don't want
z

-z

by that of

to appear in the right-hand

side of the reflection rule; also we put this in the

denominator to match the other terms in the formula.

Furthermore, we have taken account of the complex conjugation
of
z

in the law.

Geometrically it is clear that the shift in

should be of higher order than
z

-z^; otherwise we would

-44-

expect that this orthogonality-correcting procedure would

be unstable.
In addition, the other term in the reflection

rule is of order 0(h) = 0(Iz^-Zq|), and if the shift is

not of higher order, it may become the dominant term and

obscure the free-boundary condition.
term also appears to lowest order

Recall that this

in the steepest descent

equations, and that our goal is to produce a faster

convergence process.

Finally, the appearance of a tangential

derivative

(z

in the formulation should make us wary of

instability unless the term is of higher order.

These

considerations lead us to the conviction that the correction,

as formulated, should be multiplied by h
3
,

and the modified

reflection rule now reads

z.

= z +
.

-h^^^+ih.^^ ^ , i(x X +y y ) -,,3 u V 'u-^v ^v^uv + Ah ^v ^-~^0 ^_~^o

^
,

where

is a constant to be chosen experimentally.

The value of the constant

which gives good results

for all the problems considered and all mesh sizes seems to

be about 50.

However, experimentation has revealed that this

type of reflection law, coupled with an overrelaxation

technique to solve the equations, is not very sensitive to

the particular value of
A

used.

Evidently the correction

successfully achieves our goal, i.e., the enforcement of

orthogonality of
^

-^

and

v'

so that the solution is

independent of

-45-

Garabedian's basic integral formula quoted at the

beginning Df this section from reference

[3]

also provides
in theory,

a reflection rule for the axisymmetric case,

but it is much more difficult to apply the formula in a

finite-difference scheme because

the integral is not analytic.

is not harmonic and

However, an obvious modifica-

tion of the two-dimensional rule suggests itself.

Instead of
-9iP

_
'

'^v

we want
1
3ip

v'
,2

^v
l^vl

along

F;

differentiating this with respect to u, we also want

- Im
z z

vv

1 3 =-.,-2

du

,^2
y

These will be consequences of the following modification of

our two-dimensional reflection rule:
2

z,

,2 ^v ih 9 ,^v. ( _h --. 5 ~7 + -5 2 511 ^ y i = z_ +

,

'

+ Xh"

, i(x X +y y .ii.Juv-'u-'v

where we again write in a self-annihilating term which

urges orthogonality.

Reasoning in an exactly analogous

manner, we see that this rule implies the conditions

stated above. We point out here that we have just arrived at a

reflection rule by way of a heuristic reasoning process

-46-

quite different from the procedure outlined in Part

We do not know the

location of

T,

i.e. we do not know

F{z,z), and we have not invoked the analyticity of the

curve, yet we are able to use reflection to express the

boundary conditions.

This demonstrates the possibility

of applying the method to a very general class of

problems
We conclude that using the above approximations
to the reflection rules in the vena contracta problem

should insure the following three conditions:

(1^

7-72- =

''

-^

7-77
^
,

= ^

,2
(2)

-Im(z

vv V

1 8 = -^ ir 2 3u

{\p

^v

or

-Im(z

,''^v. 3 = 1 ,r (^)
55-

du

(3)

X X + y y = u-' V u V
-^

Condition

(3)

is essential for uniqueness for the

difference equations and for analyticity; as shown in

Part
1,

if it holds on the boiindary it holds everywhere,

(x
2 2 2 2
)

+y -x -y /2 must v-'v u-'u^ be a constant, which will be zero when z (w) is analytic.

and therefore its harmonic conjugate -"^

Equation

(1)

expresses the constant pressure condition

(2)

when
(1)

z (w)

is analytic, and

is just a consequence of

and analyticity.

Clearly, using a finite difference approximation

to the reflection rule requires more care than the exact

-47-

formulation.

First, we can only conclude that

(3)

holds

if our solution is insensitive, to the value of A; this

must be tested.

Second, we have no direct proof that

;

will map boundary to boundary

based on the

our argument for this is

observation that the solution to the

difference equations seems to be unique, as indicated by

the fact that the overrelaxed iterative method of solution

does converge, while we would expect it to oscillate if

there were more than one solution.

And since our equations

are certainly consistent, our solution must be the right

one, with the right boundary correspondences.

Third, we

have the problem, inherent in any reflection scheme using

four point normalization, that x +y -x -y need not be zero; ^ ' v 'v u -'u
this will be taken care of in section F.
2 2 2 2

To summarize, Garabedian's integral formula m/2 m/2 _ (z-t) (z-g(t))

-

(x,y;m)

= Im

ijj^
J

^1

^(-^,-^;l; -f
^

(^-t)(z-g(t))
)

1/2
^.^^^

^^

(i-t) (z-g(t))

provides

reflection law along the free boundary in both

In two dimensions, the

the plane and axisymmetric cases.

law can be written simply as

z

z+ - Zq " g(z_)

- g(zQ)

.2

dz

^0

but in three dimensions it is much more complicated.

-48-

Our finite difference formula approximating the rule is

given for either case by

where

we use difference approximations to the derivatives

appearing on the right-hand side.

E.

The Separation Point

When the free boundary reflection rule is applied

along the top BC of the rectangle and the y-axis rule is

used on the side AB, the corner B is naturally mapped onto

the intersection, i.e., the separation point
z

= iY-

(Fig.

5)

Since the free surface separates smoothly from the wall,

the angle between these two curves is 180, while the

corner at B is, of course, 90.

This behavior dictates

that we look more closely at the reflection rule here.

However, as we mentioned in section B, the doubling of this

angle is desirable for accuracy since ^ has an expansion

in powers of
(z-iY.)
'
.

A detailed picture of the corner B and its image

is given in Figure 7, where the points of interest w,
.
.

,w_

are indicated.

The corner B is denoted by ww-,

the neighbor-

ing mesh points are

and w^ on AB and BC , respectively,

w-^
.

and the neighboring interior point is

The reflections

of w^ and Wg through the side AB are w^ and w_

respectively.

-49-

and w_ and w. reflect w,

and w^ through BC.

The z-images

of these points, carrying corresponding subscripts, are

also shown.

The separation point is z.,

z, is

z,

lies above Zq

on the Y-axis,

below

on the free boundary T,

and z^ lies inside the flow region. The mapping

Wf., z (w)

doubles the angles at the corner

so points which lie 180 apart in the w-plane become

360 apart, or nearly coincident, in the z-plane,

For

example, w, and its reflection w- lie on opposite sides of Wq

on the v-axis, so
the y-axis.
z^

and z lie on the same side of z- on

z.

Similarly the points

and z^, and z- and z_

coalesce
The rule for reflecting through the y-axis
z^ = - i_
is exact and can be applied accurately at this corner,

yielding

z^

and z_ from

z-,

and z,.

But our approximation

T

for the rule for reflecting through

^

becomes useless;

and

which appeared in the numerator and denominator,

iJj

respectively, are both zero here, because

z (w)

on AB and

has a branch point at B

z^

A more accurate rule must be

derived for reflecting

through z^

Therefore we go back to the integral representation

in two dimensions:

^1
= ^2 - ^0
r

C'^ dz
^0

As we indicated above,

has an expansion in powers of

-50-

{z-z_)

1/2 '

More explicitly,
+ A^Cz-Zq) ^/^+. ..

= Cq + Aj^Cz-Zq)-"-/^ + A2(z-Zq)

where A^ is zero (cf

[7])

This gives

so we can write

z.

Since
z,-z is 0(h
2
)

z_

is the separation point,

3
)

z,

Thus

and the error here is 0(h

2
)

while the

left-hand side is 0(h

Also we have

but analyticity requires

A

^x

ilj

^y

on AB, so
^0 = i ^x

and the integral representation yields

^2 - ^0 = - ^^'l-^O^'^x + ^^^^

In order to use this as a reflection rule we must

find a sufficiently accurate approximation for ^

We expand
i|j

at Zq

around z^ to get 4j(x^,y^) =

ip^,

namely

-51-

since

ij;

= 1 on the y-axis

ip

vanishes and the above

estimates yield
'''3

'''O

'''x^''3~''0^

^ ^^^^

giving the approximation

4;

+ 0(h

Using this in the approximate equation for the

reflected point Z- = ^(^j^

'

^ obtain
^3-1^0
2 3

Therefore we take as our reflection law at the corner B,

'J^3-'^0

Z2 = Zq +

(Zq-z^)

{^^-^)

,2
.

Let us examine the consequences of this rule.

Since (z_-ZQ)/h and (ZQ-z^)/h are approximations for 3z/3v
at w
,

if we assume convergence of the scheme and take

*

limits as h

we get
z

= z

,2

V ^x

ip

From this we can conclude either

(1)

Im

and

^p'^

^x

= 1

or
(2)

z^ =

But

x=OonAB, SOX V =0
'

and either conclusion leads to

=0

-52-

With this in mind we divide by h 2 /2 and take the

limit.
^

We get ^

vv

= (z
'

ili

^x

Now we can conclude Re

z 2
<iJ

vv

=0

(which is no surprise) >^

'

and either Im
z

vv etc., etc., until we eventually get a non-zero derivative

/

In the former case, we have

=0

=0 also or X =1. vv 3 h /3' r and we divide by J

2
ii

'V'

we will finally conclude at which point '^

-^

^x

=1.

Therefore we see that the consequences of this

reflection rule are

z

and

2
i|^

^x

=1.

The first condition is clearly consistent.

Since

^p

is the normal derivative of ^ in the z-plane at

the separation point, the second equation is the constant

Speed condition, and thus it is consistent too. Summarizing, we have shown that in the two-dimensional
case the reflection rule for obtaining z^ = zCw-)

from z^

should not be the same rule as was used along BC; instead
it

should be calculated from

z, = z +
2

(z^-zj 1

-)'^ { ^x^-x-.

In three dimensions, this is changed to

Z2 = Zq +
'^B'^O (^)2
2

(Zq-z,)

1^
y

Entirely analogous calculations show that this law implies

y which is the constant speed condition at separation.

-T

'^x

= 1

'

-53-

In closing this section we add a note about the

implementation of this theory.

One desired effect of this

reflection rule is to cause the coefficient of (z.-z^),

which approximator, the speed of the flow at

unity as h goes to zero.

to approach

However, since the speed is

already forced to assume this value along the free boundary

by the other reflection rule, one might assume that this

coefficient may well be replaced by

in the formula,

merely allowing continuity to enforce the constant speed

condition.
The results of numerical experiments confirm

this reasoning, and we found that comparable accuracy could

be achieved by the simple rule

Z2 = Zq + Zq - z^

in both

the planar and axisymmetric cases

F.

The Line of Truncation

The values of x, y and

i^j

on the vertical line CD

truncating the infinite rectangle in the w-plane will approach

their limiting values at infinity as the length of the

rectangle, AD, hereafter called U

As an approximation we set
y = v
,

max ,

is increased

(see Fig. 5)

X = x.

= constant,

and
=

ij,

-54-

'

on CD, and solve the resulting set of difference equations

for several values of U U

max ,

then extrapolate our data for

'^

max

-*

To do this we must decide how to select x,

Recall that no provision has yet been made for

2 2 2 2 driving ^ the constant H = x +y ~x -y V 'v u -'u

to zero.

Since

x, 1

is

the only undetermined parameter left in the problem, it is

here that we enforce this condition.

Experimentation on

several models has indicated that the solution to the

difference equations is unique for any fixed

that H is a function of x,
.

x,

and further

Therefore it is through our

.

choice of

x,

that we make H =
is too large,

If x-

the stretching in the u-direction

,

of the map

z (w)

is greater than in the v-direction

so H H will

will be negative; on the other hand if x^ is too small

be positive.

We can use this knowledge to calculate

x,

while

we are iteratively solving the difference equations.

calculate H at each iteration, then shift
6

We could

x,

by an amount

X,

= y H

where

y is

a positive constant.
[1]

This is essentially the

scheme in

The scheme which we shall use is a slight variation

of this which is faster and more stable.

More specifically,

H is calculated by finite differences at every point in the

w-rectangle to the right of u =

1,

and an average is formed.

Although H should be constant once the solution to the

-55-

difference equations has been achieved, it need not be

constant during the iterations, and this averaging makes
the procedure more stable.

Since we shall exclude from

our average the part of the rectangle containing the pole

and the branch point, the process is stabilized even more.

We calculate 6x, by the above equation; however,

we not only shift x

by this amount, but all the x-coordinates

1

of the mesh points to the right of u =

are shifted by a

proportional amount.

This not only stabilizes the scheme,

but also speeds up convergence.

Finally, rather than apply this procedure at each
iteration, we wait an appropriate number of iterations after

shifting

x,

for the effect to settle, so that the next

calculation of H is more meaningful and realistic.

also increases stability.

This

Experimentation indicates that

the optimal number of iterations performed between shifts

of x^ should be about equal to the number of mesh points

along ^ the u-axis

, '

or in other words U

/h. max'
y

We can estimate the optimal value of the constant

by the following argument.
a
(4a) xl

If a 4x1 rectangle is mapped onto

rectangle using our reflection laws,

the mapping

is given by

X = au

y = V

This gives
H = x

2^2 +y

-'v

2 2 -y -X u -^u

=l-a
T

and
x,

= 4a

-56-

The mapping is conformal only when a = 1, which means x,

must be shifted by
6x,

- 4a

Thus

ought to be the ratio

5^1

^ -

"T" 1.

4 _ 4- 4a _ J - T+E : 1 -a

giving

y =

when

a =

We anticipate that when our w-rectangle in the vena contracta problem is 4x1, the optimal value of
in the shift formula should be aro\ind 2.
y

Actual

experimentation

shov/s

that the best choice is

y

= 2.8

G.

The Singularities at the Origin

As we mentioned in Section B

the corner A of the

rectangle in the w-plane, located at the origin, is mapped

to the flow region at infinity in the second quadrant of
the z-plane (Fig.
5)
,
.

This requires a pole in the

so that near w =
^

transformation

z (w)

z(w)

where R is a negative number.

This behavior naturally produces an unbounded

truncation error in Laplace's difference equation, so it

must be compensated for in some way before the problem can

-57-

be handled on the computer. to that in [1]

.

The method we use is similar

In the triangular region of the rectangle

determined by the base AD, the side AB, and the line
running from
vr

= i/2 to

w = 1/2

(Fig.

8)

we subtract off

the pole and solve for the function

= z(w)
-

?(w)

^ w

This function is still analytic; thus it too can be

approximated as a solution to Laplace's difference equation.

Furthermore
region.
?;

(w)

has no singularity in the triangular

Therefore the truncation error is bounded and

a machine solution is feasible.

Clearly, knowing
to knowing
z (w)
.

^ (w)

and the residue R is equivalent

It remains for us to specify the boundary

? (w)
,

conditions to go with the difference equation for

to calculate R.

and

We will discuss the boundary conditions first,

? (w)

There is no problem in specifying the values of

along the hypotenuse of the triangle, because we have values

of
z (w)

at these points; we simply subtract the values of R/w,

The details of the interfacing are quite straightforward,

and we omit them here (cf.

[1])
C

The boundary conditions for

along the two

sides

of the triangle enclosing the corner A are derived from the

reflection laws for

the law for
z

z.

Since AB maps into the y-axis,

is
z^ + z_ = +
z (w_)

z (w_^)

-5 8-

where

w=-h
the

iv,

=h+iv
size and v denoting

with h, as usual, denoting the

rr.esh

ordinate of the mesh point.

Substituting

r,

+ R/w for z and noting that R is

c

real, we see that the reflection rule for

C+ + ?_ =

is given by

i.e., it is identical with the rule for z.

Similar considerations
reflecting
r,

shov7

that the law for

along the bottom of the triangle is

?+ - C_ =
.

The procedure for determining

? (w)

is now completely

specified except for the calculation of the residue R.

We turn our attention to this problem.

First we point out that if the function

to be analytic, the residue R is determined.

z (w)

is

We have
1

completely specified the flow region by setting Y^ =

(cf
,

Section

A)

and we have fixed the images of three

boundary points on the semi-infinite rectangle in the w-plane

(the point at infinity and the two corners)
,

so there is

no more freedom left in the choice of the confoonnal map.

Numerical experiments were performed, solving the

system of equations described above v^ith fixed preassigned
values for R to see the effect of this parameter.
It was

found that the iterates converged to solutions wherein the

-59-

free boundary detached from the fixed wall at some non-zero

angle, rather than detaching smoothly as the theory predicts.

Furthermore, in the neighborhood of the separation point

the Cauchy-Riemann equatioTiS were ivolated.

However, we

noticed that this angle at detachment was a monotonic function

of the parameter R, so that any desired angle could be

achieved by the proper choice of the value of R.

This immediately suggests an iterative procedure
for calculating the residue; one starts with an initial guess and then adds to this value a correction depending on the

angle formed by the free surface and the wall at detachment,

until this angle is zero.
The iteration of the value of R

can proceed simultaneously with the iteration of the

difference equations.
this schem.e.

There are many ways of implementing

We will report the details of the method

as judged by its success in

Vvhich seems most accurate,

the two-dimensional case.

We refer again to Figure

7,

where the mesh points

,

along the top edge of the rectangle are called w

w^-

Wq,

and Wq and their z-images are similarly designated.

The

smooth detachment property at the separation point is

demonstrated analytically in the power series of

the top edge, expanded at w = w
.

v/

along

Denoting the abscissa

by u (the ordinate is

i)

we have

-60-

x(u+i)

^0

-^

y(u+i) =

at separation:

The constant

is chosen experimentally to achieve optimal

convergence, but once the system has been

solved, the

solution should be insensitive to changes in n, since x

must be zero.
The best formula for calculating x
in terms of

differences involving neighboring mesh points is

9x
=

(3xg -

I Xg

^ Xg)/h + 0{h^)

This is the exact formula for the derivative of a cubic

polynomial, so its use is most appropriate here, where

we seek the dependence in the form

x near separation.

'^

We can summarize our procedure for calculating the

residue R as follows.

It was observed that the calculated

angle between the free boiondary and the fixed wall is a

monotonic function of the value used for

R.

Since the

theory states that this angle must be zero, this property

is used to find the correct value of R.

An accurate measure

of the deviation of this angle from zero is provided by

the derivative x

as approximated through the above

difference formula, and this value is used in evaluating

the corrective tern

-62-

5R =

5^

in an iterative scheme for calculating R.

This iterative scheme proceeds along V7ith the

iteration of the difference equations for

and

ii

but rather than correcting R at each cycle, we allow

the effect of a correction to settle down before correcting

again.

In fact, the iteration of R is done simultaneously

as described in Section C.

with the shifting of x

We add that the above process works satisfactorily

for the vena contracta, and that when the iterations have

converged the Cauchy-Riemann equations hold at w

cubic approximation for x

The

gives significantly higher

The optimal value for
r\

accuracy at every mesh size.

seems to be about 30

We conclude this section with a brief discussion

of the equation for
i|^

near the origin.

rP

Since

= 1

on AB and
li;

on AD, this discontinuity will affect the accuracy of the

difference equation unless it is properly handled.

The treatment here proceeds exactly as in
[1]

In the triangular region we subtract off a solution ^,,

which incorporates the behavior near the origin, so that

the remainder has continuous boundary values
,

In two

dimensions.

-63-

ill,

2 arccos

,r

TiT+vV^

t-t-tt'

and in three dimensions,

I

''^1

U
,,

u
2

1/2'

In tJiree dimensions the truncation error in the

difference equation for

becomes unbounded at the origin,

and this, too, is handled by a modification inside the

triangular region.
to
z'

The variables

and

ip

are transformed

and

ip

'

via the Kelvin transformation

z'

= 1/z =
i(;/lz|

^'

Then the equation for ^'{x',y') has the same form as

that for
'|'(x,y),

and the truncation error is again

[1]

bounded,

For details see

H.

The Interior Points

?.s

we said before, the equation

Ax =

is

approximated by Laplace's difference equation

X.
.

1,3

T
TT

(X.

+ X.

1-1,

i,D-l

+ X.,,

1+1,3

+ X.

ifD+1

.,,)

written at every interior mesh point.

is

A similar equation

written for
Ay =

-64-

and, in two dimensions, for

Alp

In three dimensions, the equation for the stream

function, which we write in the conformally invariant form

can be approximated to the same degree of accuracy by

(cf.

[ID
=

i-l/J J--1/3
a.

1,3-1^1,3-1
i-lfD
,

1+1,3^1+1,]

i,j+l^i,3+l

,|,.

'-'

.+a.

i/D-1

,+a.., .+a. .,, i+l/D 1,3+1

where
^'^

^k,+yi,j

This equation is written at the interior points in the

axisymmetric case.

After the Kelvin inversion mentioned

in the last section, the equation retains the same form

for the transformed values.

These finite difference

equations are handled by the standard procedure of

successive overrelaxation, which seems to be the best

available vjhen nonlinearities are involved.

-65-

I.

THE FORTRAN PROGRAM

In Appendix II we list the Fortran program used

to calculate the vena contract as described above.

A typical

run involving 900 iterations on a grid of 11,000 mesh points

required 40 minutes of machine time.

However, due to the

operating system on the CDC 6600, the same result could

have beer achieved in about a third of this time if the

unknowns x = X, y = Y, and

i|j

= P had been stored in

single

instead of double arrays.

We also mention that when a

similar large run was tried on the IBM 360/65 at the

University cf South Florida it failed because of the smaller

word size.
Careful study of the code will enable the reader
to identify
tJie

computations described in the text, but

it will be helpful if we define the input paraur.eters

M, N

are the number of mesh points in the v and u

directions respectively

NOHALF
IPOLE

IS the number of times the mesh is to be halved

is the location, on the v axis, of the interface

defining the triangle at the origin (Section

ITPSIO
is the number of iterations

G)

performed on x and y,

per iteration on ^

ITERMAX is the maximum number of iterations allowed

EPS
CON
is m, where m+2 is the number of dimensions

is

of Section A, Part

-66-

RELAXZ, RELAXSI

are the parameters c in the relaxation

^

factors r =

2 ^-^

in the computation of

and

ii

respectively
CONVERG is the number with which the residuals are compared
to determine convergence

RESCON
SENS

is n of Section A, Part

is y of Section A, Part 4.

-67-

Part 4.

RESULTS
The program described above was run on New York

University's CDC 6600 computer for many test cases.

Here we shall describe the results with regard to the

parameters, accuracy, and convergence.

We conclude with

our new estimates of the contraction coefficient.

Throughout Part

we shall compare our results

[1]
,

with those of Bloch in reference

because the m.easure

of our success is the degree to which the reflection scheme

improves on the efficiency of the method of steepest

descent.
In Section C on rates of convergence the

superiority of the reflection technique is clearly demonstrated

by computer runs requiring about one-tenth the number of

iterations needed for analogous cases employing steepest

descent.

-68-

'

A.

Dependence on Experimental Parameters

There are five parameters in our vena contracta

model which are chosen experimentally to achieve certain

goals.
(1)

These are:
the coefficient of ^^^^'''^u^v ^^ ^^^ free

X,

boundary reflection rule, which enforces orthogonality,

(2)

y,

2 2 the coefficient of x +y -x^-y^ in the equation

for shifting x,

the truncated end of the jet.

2

The goal here is to nullify ^v+yv~Xu~yu

(3)

n,

the coefficient of

at

separation, which occurs

in the equation for correcting the residue R in

order to achieve smooth detachment.

(4)

the location of the hypotenuse of the triangle

bounding the region containing the origin of the

w-plane; in this region singular terras are

subtracted off from

z (w)

and

ip

(w)

In the text

we choose this line to run from w = i/2 to w = 1/2, but

(5)

other choices are peinnissible.

the length of the rectangle in the w-plane

max ,

This rectangle m.ust be long enough to produce

accurate approximations to the theoretical results

of the infinite jet.
It is clear that the rate of convergence of the

iterates of the overrelaxed system of equations will be

affected by the choice of A,

y,

and n

but the final solution

-69-

must be insensitive to them if they are to be effective

in achieving their goals;

i.e., the terms which they

To get fastest convergence,
A

multiply must go to zero.

we begin the iterations with the values

n = 40.

= 50,

= 2.8,

Once convergence is achieved,, as indicated by all

residuals being less than some fixed number, usually 10"

we double and triple each of the numbers X,

and

n to

test the insensitivity of the solution to them.

In the

cases tested, this perturbation required only a few more

iterations (less than 10% of the number of iterations

needed for initial convergence) to reproduce the convergence

criteria, and the values of the contraction coefficient were

changed by no more than one unit in the fourth decimal place.

The location of the interface between the triangular

subregion at the origin and the rest of the rectangle should

have very little effect on the data for small mesh size,
since the term being subtracted is an exact solution of the

appropriate differential equation.

Testing this witli

different choices of the interfacing line produced no change

in the fourth decimal of the contraction coefficient nor in the rate of convergence of the iterates

max can be estimated by performing the computations for several

values of U^ ^ and extrapolating for U
*

The effect of truncating at u = U ^ the rectangle ^

In fact,

Bloch

[1]

has argued that the extrapolated contraction

-70-

coefficient, C
C (U
)

differs from the approximate va3.ue

according to the asymptotic formula

C = C
(U

max'

+ A e"'^ "n^ax
.

where A and

are constants

We shall examine this

phenomenon carefully in Section D, but first we turn

our attention to the accuracy and convergence of the method
on a rectangle of fixed dimensions,

B.

Accuracy of the Solutions

In this section we shall verify that as the mesh size

h is decreased, the contraction coef f icierit, which is

representative of the values of the solutions on the mesh

points, approaches a limit with an error varying as the

square of the mesh size.

We anticipate this kind of behavior because

throughout the program we use approximate formulas whose error

are no bigger than
2 0(h).

For example, our reflection

v/e

rules are all at least this accurate, and

write centered

differences or high order extrapolated differences for

derivatives.
The data presented herein supports our
)

contention that the accuracy is 0(h 2

We tabulate the values of the contraction coefficient

calculated at the various mesh sizes in two dimensions, with

Vx

=4

1/3.

-71-

C.

Rate of Convergence of the Iterations

It was expected that the iterates of the reflection

schene would converge much faster than those of the method

of steepest descent.

The computer

rvins

confirm this, and

we present here a comparison of two of our runs v/hich

correspond closely with the runs quoted in
[1]

In the two-dimensional case we chose a rectangle

with
U

max

=4

1/3, ' '

h = 1/48 '

It involves 10,192 interior points and 516 boundary points,

totaling 10,70 8 mesh points.

The computer
;

V7as

run until

all residuals were less than 10~

the iteration was started

with the exact solution.

the run was
5 88.

The number of cycles required for

For comparison, Bloch quotes a two-

dimensional run starting with the same data and cutting off
at 10
6
;

it employs U

max

, '

h = 1/48 '

and 14,161 mesh points.

The run requires 6,260 iterations.

For the axisymmetric case the reflection program was

run

with U

max

=4

the 1/3 again, so that V7hen h = 1/48, w w ^ /

I

number cf mesh points was 10,70 8 as before.

However, to

initialize the iterations the following procedure was used.

The exact solution in two dimensions was computed to start
the procediare, with a mesh of h = 1/6;

after residuals v?ere

reduced to 10

the mesh was halved and the data interpolated

-73-

for the new mesh points.

This was repeated, but for

h = 1/48 the iterations proceeded until the residuals

were reduced to 10~

The nuriber of iterations required

at each mesh size is presented in the following table:

Niimber

D,

Estimation of the Contraction Coefficients

As mentioned in Section A, the calculation for

the contraction coefficient should be corrected for the

truncation of the rectangle by a formula (cf.

C
= C
(U

[1])

max'

+ A e

max

This rule should be valid at any mesh size, with A and

independent of h, at least when h is small enough.

We tested this behavior by computing C values of U
for different

max are as follows

u

The results for the mesh size h = 1/24

max

The values for

differ only insignificantly from

The coefficient A is of

those of [1], as we should expect.

the opposite sign, however.

This is probably due to the

fact that Bloch uses different boundary conditions for y and

^

on the truncating

line CD.

To estimate the contraction coefficients, we start

with our most accurate values, computed for

U

max

=4

1/3, r /

h = 1/48 /

and residuals less than 10~

-

We first extrapolate these values for h

'

by the h
--

formula of section B, then extrapolate for U ^

the above.

max This leads to the final estimates

C
=

>

by
-^

.61106

m = m =
1

C^ = .59142

where at least the fourth decimal is significant.

For comparison, Bloch
[1]

quotes

C^ = .61100 + .00002
C = + .00004 .59135
, '

and our computations confirm these values

In closing, we add the following observation.

The coefficient in the

planar vena contracta is calculated

from the exact solution to be

=

.611015...

so our estimate is high by .000045.

Therefore a prudent
in the axisymmetric case

guess for a five place value of C

is obtained by subtracting this same correction,

yielding

^c = .591375.

-76-

BIBLIOGRAPHY
1.

Bloch, E., A Finite Difference Method for the Solution

of Free Boundary Problems
^

AEC Computing and Applied

Mathematics Center, Courant Inst. Math. Sci., New York

Univ., NYO-1480-116, 1969.
2. Garabedian, P.
R.
,

Partial Differential Equations ,

Wiley, New York, 1964.

3.

Garabedian, P. R., Calculation of Axially Symmetric

Jets and Cavities
,

Pacific J. Math.

6,

1956.

4.

Garabedian, P. R., Estimation of the Relaxation Factor

for Small Mesh Size , Math. Tables Aids Comp. 10, 1956.

5.

Garabedian, P. R. and Spencer, D.

C,

Extremal Methods

in Cavitational Flow , J. Rat. Mech. Anal. 1, 1952.

6.

Giese, J. H., On the Truncation Error in a Numerical

Solution of the Neumann Problem for a Rectangle ,

J. Math. Phys.
7.

37,

1958.
,

Gilbarg, D. Jets and Cavities

Encycl. of Physics 9,

Springer-Verlag, 19 60.
8

Lewy , H

. ,

On the Reflection Laws of Second order Differ-

ential Equations in Two Independent Variables ,

Bull. Amer. Math. Soc. 65, 1959.
9. Nehari,

Z.,

Conformal Mapping , McGraw-Hill, New York, 1952.

10. Southwell, R. and Vaisey, G., Fluid Motions Characterized

by "Free" Streamlines
Ser. A 240, 1948.
11. Young, D. M.,

Phil. Trans. Roy. Soc. London,

Iterative Methods for Solving Partial

Differential Equations of Elliptic Type, Trans. AMS76, 1954

-77-

WALL

FREE SURFACE

FIGURE
(Plotted

from computed data for the

axisymmetric case)

-78-

U
.Wa

Figure 2

-79-

Figure 3

80-

Figure 4

-81-

 z+

Zo

z-

Figure 6

-85-

W, o

w.

W7

(B)

w.

Wn

w.

w,

w.

z,

<.

Z3

Zo

Z4

z,

Figure 8

-85-

LT

APPENDIX

1.

^*^^

PROGRAM REFLCTR (I NPUT OUTPUT INTEGER PPS DATA COMMON X(405,5?).YM05,52).P(405,52).M. CWY(40,40) ,72(40.40) JST , B6S, H, CN, POLE, PPS, II. I?. ITER,ERX, IX. JX.FRY. Y JY PRS! . IS ITPSI0.8BNS CITPSI,SR.GABY4.SRP.PABY4,C0N. I, J. COMPLEX 7 read l.m.n.nohalt, pole* itpsio, terh ax. ep^. ton rblax? , pel axsi. cconverg.rescon.sens READ l.L.JUMP CALL SECOND(T) PRINT 141, PRINT 70. rORMAT (IH HO,* L ) 70 Ml./rLOAT(M-l) IlaN*l I2=M*2 UMAXN*H 1 FORMAT (7I10/7E10.4) CALL INIT 20 PRINT 2,M,N,IP0LF. ITPSIO. I TERK AX FPS CON. 9il AX? REL AXST , rONVERa, CRESCON.UMAX.SENS.NOHALF ITERMAX , N ITPilO M IPOLi 2 FORMAT (#0 RELAX PS RELAX Z /IH ,6I10/*0 CCNTRCL EPS C */lH .^Fl9.4/ CONVFRQF RESCCN CI*,* SENS NOMA F */ UMAX C*0 CIH .?E15.4. HO) IWAVE=40
, .
I I . .

ITER-0 JERRMAX100. $JPRNT1 IPRNT50D TEST X in. ITPS10 CALL PRNTR(NOHALF, SHIFT)

SHI FT. 11 111111

.

SRs2./(l.*M*RELAX7) QABY4a.29*SR SR l.-SR PABY432./(1 .*H*RELAXSI SRPbI .-PABY4 IF (EPS.PO.O) PABY4.25*PaBY4 ITER ITFR*1 CALL REFLEK
)

CALL INTRR I3N-i DO 80 J3,M Y( 11, J)Y( 13, J) 80 P( II, J)aP( 13. J) ITER-IPRNT) .LT.O) GO TO IF IF (TEST.LT.ERRMAX) GO TO 60
(

-86-

PAGE

1
0

99

CALL PRNTR(NOHALr, SHIFT) PRINT 61 ERROR HAS GROWN. *} FORMAT (* CALL EXIT ERRMAX = TEST IPRNTIPRNT*500 CALL PRNTR(0. SHIFT) ITER.LT. ITERMAX) QO TO 99 IF CALL PRNTR(0, SHIFT) CALL EXIT TESTAMAX1(ERX.ERY,ERSI) IF (TEST.LT.CONVERQ) GO TO 11 ITER.LT. IWAVE) GO TO 4 IF CALL COMPMOD(SHIFT.RESCON, WAVE, L)
(

QO TO 4
7

IF ({ ITER-JPRNT) .LT.O) GO TO 9 JPRNT JPRNT JUMP CALL PRNTRd. SHIFT)

00 TO 9 11
IF (ABS(SHIFT) .LT.CONVERQ) GO TC 40 CALL COMPMOD{SHIFT.RESC0N, 1WAVE,L)

QO TO 4

CALL PRNTR(NOHALF. SHIFT) PRINT 31, SHIFT FORMAT (#0 LAST SHIFT WAS SI CALL SECOND(T) PRINT 141, IF {N0HALF)12,12 .13 13 CALL HALFER
40

,E10.3)

NOHALFsNOHALF-l PRINT 14 14 FORMAT (*0 MESH SIZE HAS BEEN HALVED.*) CALL SECOND(T) PRINT 141, TIME IS ,F10.3) 141 FORMAT ( READ 50 CONVERQ. ITERMAX, RELAXZ, SENS. RES'^QN FORMAT ElO 4 110 SEIO 4 50 LL*L
, (
.

QO TO 20 12 CALL EXIT

END

SUBROUTINE PRNTR JAX, SH IFT)

(

DATA INTEGER EPS

-87-

paQE

COMMON X(405.5?).Y(405.52).P(405.52).H. CWY(40,40) ,72(4n,40). CN, IPOLE.FPS, Il.I9.ITER,ERX. IX, JX.FRY, lY. jy.PRST, ISl JST , PES. M, C1TPSI.SR.GABY4.SRP.PABY4.C0N.I,J, itpsto.ssns PRINT 1, ITER. ERX. IX, JX.ERY. lY.JY.ERSl. ISI. JST, CSHirT,RES V*/ X*/inX S50 10 15. 15. 1 FORMAT (0*# I9.E90.10, 15, 15,* PSI*/ C10X,P20.10, 15, 15.* SHirT*/lH . B20 10 . ?flX. BES* C10X.E20.10,10X. CC(1 ./Y(2.M*1) )**(1*EPS) PRINT 30. CC F11.7) CONTRACTION COEFFICIENT FORMAT (0 SO PRINT 300, X(N*1.?) XCEND)" ,F10,6) 300 FORMAT ( RETURN END
.
.

SUBROUTINE REFLEK DATA COMMON X(405,5?).Y{405,52).P(4 05,52).M.

CWY(4O,40),Z2(4O.4O).
POLE, EPS, 1. I?. ITER, ERX, X JX ERY lY. JY.PRST , TSI . JST .RES, CITPSI.SR.GABY4.SRP,PABY4,C0N, I.J. ITPilO.SBNS
ON,
I
1

H.

INTEOER EPS COMPLEX 7.RHS IPMIN2= IPOLE DO 1 I-IPMIN2.N X( I,1>hX( 1.3) Y( I,1)B-Y( 1,3) jM*l DO 2 I" IPMIN?,
X(l,

Yd,

)"-X{3, 1) 1)Y(3, I)
1

IPMIN2IPMIN2-1 JIPMIN2-1 DO 3 I=J. IPMIN? Zb RgS/(CMPLX(FLOAT( I-2)l. )-) X( 1,1) X( 13) -REAL(Z) Y(I,1)3-Y( 1,3) A1MAQ{Z) 7s RES/(CMPLX(1. .FL0AT(I-2>)*H) * REAL(Z) --XCS. X(l, Yd, I):Y(3, I)-AIMAO(Z) JJ-1
1
)

DO

1=2,
1

X( I,l)aX(

.3)

-88-

PAQE

V(

Ijl)n-V(I,3)

Yd, I)Y(3.I)
IF

(EPS) 5.5.6 5 DO 51 I3.N C-( .5*P( I.J-2) -2P(I#J-1) 1.5P(I.J))**9 )-X(I.Jin D (X( I*1,J)-X( I-1.J))*(X( I,J )*(Y( 1*1. J)-Y( I-1#J) I, >Y(!.J-in (Y( J C*

Db.5*H*C0N*D
C

.2B*(( ,5*P(I*1. J-2Ul.5*P(I*l,J)-2.*P{

!-l.jn**2)
))

1*1 ,J-1

))*2

C-( .5*P( I-l, J-2>-2.*P(I-l#J-l)*l.B*P(

RHS =CMPLX(C,D) 2 CMPLX(X( I. J),

Yd, J))
.7>

ZC0NJQ(7)*RHS/(CMPLX(X(I. J-1),Y(I,J-1 X(I,J*l>cRBAL(Z) 51 Y( I,J*1)-AIMAQ(7) R(P(3,M)-P(2,J)>/(X(3,M>-X(?,J)) 52 ZbCMPLX(X(2.M>-X(2,J),Y(2,M)-Y(2, J)) 2CMPLX(X(2,J).-Y(2,J)> -Z*R**2 X(2,J*1) - R6AL(Z) Y(2,J*1>-AIMAQ(7) X(3,J*1)X(1,M) Y(3, J*l>Y(l,M) RETURN
6

)-X(I,J-l)) C*<Y<1*1,J)-Y(I-1.J))(Y(I.J )-Y(I,J-l)) PV.B*P(I,J-2)-2.*P(l, J-1)*1.5*P(I,J> YU,5*(Y(I*1.J)-Y(!-1,J)) PUV.25*(P( 1*1, J-2>-P{ I-l, J-2) >-(P( 1*1 , Jn-P( T-1, J-l)> C*.75*{P(I*1,J)-P(I-1,J)) YIN1./Y(I,J) )**9 C-{PV*YIN DONaCON*M 90 D.5*D0N*D(PIJV*PV*YU*YIN)PV*YIN*2 RH8CMPLX(CD) 91 2CMPLX(X(I,J),Y(I,J)) 2C0NJQ(Z) * RHS/{CMPLX(X(I,J-1).Y(I.J11) -7)
X( I, J*l) n REAL(7) 61 Y( I,J4>1)BAIMAQ(7)

DO 61 183, D=(X{ I*l,J)-X(I-l.J)>*(X(I,J

R(P(3.M)-P(2.J))/((X(3,M)-X(2,J))Yep,J)> ZCMPLX(X(2,M)-X(2,J),Y(2,M)-Y(2.J)) Z*CMPLX(X(2,J).-Y(2,J)) -Z*R**2

X(2,J*1) REAL(7) Y(2, J*1>-AIMAQ(7) X(3,J*1)-X(1.M) Y<3,J*1>Y(1,M)

-89-

PAGE

RETURN END

SUBROUTINE INIT DATA COMMON X(405,52).YM05.92).P(4 05.52).M.

CWY(40,40),72(40.40), CN, POLE, EPS, II, I?, ITER,ERX, IX, JX,eRy,IY,jy.PRST. !8IiJST,B6S.H, C1TPSI,SR.GABY4.SRP,PABY4,C0N. I.J. ITPSlO.ifNfi
I

INTEGER EPS COMPLEX W. ZETA, E* 62. 8. M1M*1

L,

DO 1 1-2, DO 1 J2.M1 Us CMPLX(FL0AT(I-?),rL0AT(J-2))*H

2

IF {I*J-4) 2,2.3 X(2,2)30. Y(2,2)aO. P(2.2)>0.

00 TO
3

ZETA0l.57O79633*U ZETA .5(CEXP(ZFTA) - CexP(-ZBTA)) ZETA* .63661998 *CLOQ(ZETA) P(I,J)* AIMA0(7ETA)

E

-1.57079633*ZPTA

E>CEXP(E) E2"EE*1. S-CS0RT(E2) GAIMAG(S} IF (Q.LT.O.) SCONJQ(S) LCL0G(1.*S) OaAIMAG(L) IF (Q.LT.O.) LC0NJG(L) Zs ZETA*. 63661 998* (E-S*L) X(I,J)-REAL(Z) Y(I.J)=-A1MAG(Z)*2. CONTINUE N1bN*1 X(N1.2)3X(N,2)*H Y(Nl,2)n. DO 4 J=3.M1 X(Nl.J)aX(Nl,2) Y(Nl.J) (J-2)*H IF (EPS) 5,5,6 PCNl.J) Y(N1,J)

-90-

U
PAGE

QO TO
6 4

P(N1,J) Y(N1,J)**2 CONTINUE P(N1.2)0. IDIP0LE-4 RES-(l.-EPS)*.8lQ97*EPS*(-.6055) IF (EPS.EQ.O) 00 TO 8

DO 10 I?. DO 10 J2,M1
P(

10
8
7

I,J).B*P( I.J)
II3, ID

DO 7

CALL FLIPZ(-1. 11) CALL FLIPSK-l. II) RETURN BND

SUBROUTINE FLIPSKK.ID)
DATA COMMON X(405.52).Y(405*92)P(40592).H. CUY(40. 40). 72(40, 40) CN, IPOLE.EPS, Il.I9.ITER,ERX.lX,JX,FRY.TV.Jy.PRST,TSI*JST,RES.H.

CITPSI.SR.QABY4.SRP.PABY4.CON.I,J.ITPST0.SeNS
INTEGER EPS Ie2.ID JID-I*2
DO 3
2
1 3

WS0RT(FL0AT((l-?)**2 *(J-2)**2)) IF (EPS) 21.2 P(I,J) P(I,J) .5*K*(l.-FL0AT(!-2)/W

00 TO 3

P(I,J)sP( I,J)* K*2*(l.-.3l830989*AC0S(FL1ATn.2)/W)) CONTINUE RETURN END

SUBROUTINE FLIPZ(K. ID) DATA COMMON XM05,5?).y(405.52),P(4 05.52).M, CWY(40, 40), 72(40,40), CN, POLE, EPS, U, I?, ITER,ERX, X, JX FRY I Y. JY PRST CITPSI,SR.QABY4.SRP.PABY4,C0N, I,j. ITPSTO.SCNS INTEGER EPS COMPLEX 7,W
I I ,
, ,

ISI

JST

RES, M,

-91-

paQE

DO

1=2. ID
( (

J1D-I*2 WCMPLX( I-2)*H, J-2)*H) 7 = CMPLX(X( I# J) YM. J) X( I, J)=RFAL(7) y{ I. J)s AIMAG(7) RETURN
. )

KRES/W

END

SUBR OUTIN E COMPMnD(SHIFT,RESCCN# IWAVS.LI COMP LEX 7 COMM ON X( 405,52).Y(405,92)iP(405.5?).M. CUY(4 0,40) .72(40.40). CN.IP OLE,P PS. 1. 12. ITER.ERX. X JX f RY I Y. JY PRSf CITPS I.SR, GABY4,SRP.PABY4,C0N.I,w. TPT . SBNs IWAV E=IWA VF*N SOD 0.
1

ISI

JSf .RBS, H,

90

J3 .H 1*1. XUX J)-X( 1-1. J) XV"X (I.J* 1)-X( I. J-1) YUY 1*1, J)-Y( I-l, J) YVY (I.J* 1)-Y( I, J-1) SODb SOD*( XU*YV)*(XU-YV)*(YU*XV)*(YU-XV) SOD SOD/( 4.*H**2) SODb SOD/F LOAT( (M-?)*(N-L*1))
(
(

DO 2 DO 2

IL .N

11

MleM 1 SHir T-SP NS*SOD IF ABS(S HIFT) .GT.H) SHIFT"SIGN(1.,SHIFT>H FACT 0R = 1. SHIFT/ X( N*i,2)
(

100

DO 1 00 I- L.U DO 1 00 J 2. Ml X(I. J)=FA CTOR*X( I.J) JM1

SLOP Ex-FL0AT(M-l)*(3.*X(S. J)-1.5*X(4. J)*,^3:SS*X(9.J)) Da-R ESC0N*SL0PE ABS(D) .GT. .09)) DbSIGN(1..D)*.0S IF RESb RES-D DD/ M I?.N DO 3 INTF CsMAX0(2. IPOLP-I-1) JsINTFC.Ml DO 3 ZsD/ CMPLX{ FLOAT I -2 FLOAT J-2 ) X(I. J) = X( I, J)-REAL(Z) Yd. J)=Y( I, J)-AIMAQ(Z)
(

-92-

PAGE

30

CONTINUE RETURN
END

SUBROUTINE HALFER INTEGER EPS COMMON X(40 5,5?).Y(405.52).P(405.52).M, CWY{40, 40), 72(40,40), CN, POLE, EPS, 11,12, ITER,ERX,IX,JX,PRY. TY.JY.PRST, SBNR CITPSI,SR,QABY4.SRP.PABY4,C0N. TPSI J II-lPOLE-4
I
I

IS

JST

RES. H,

DO 22 22

IIIa3, II CALL FLiPZd, III) CALL FLIPSKI.III) M1M*1 DO 1 II2,I1 IsN*3-II DO 1 JJ2,M1 JsM*3-JJ

21

X(2*I-2,2*J-2)X(I,J) Y(2*I-2,2*J-2)Yri,J) P(2*I-2,2*J-2)P( I, J) P(2,2)s.5*(2-EPS) N2bN*N M22*M-1 DO 2 I2,N2.2 X(I,3)=.75*X(I.4)*.375*X{I,2>-.125*X(I,6) Y(I,3)s.75*Y(I.4)*.375*Y(I,2)-.125*YM,6) P( I,3)s.75*P( I.4)*.375*P(I,2>-.125*P(I.6> DO 21 I2,N2,2 DO 21 J5,M2,2 X( I, J)s.75*X( I,J-1)*.375*X( I,u*l>-.J99*X( Y( I, J) 8,75* Y( I, J-1)*,375*Y( I v>l )- IS^^Y P{ I, J)s.75*P( I. J1)*.375*P( I, J*1)-.125*P( P(2,2)bO. MM2 NN2-1 M1M*1 DO 3 Ib3.N,2
.
.

T*.l-3) t ..1.3) T..1-3)

DO 3 Js2,Ml X(I, J)=.5*(X( I-1.J)*X( 1*1, J))

3

Y(I,J)=.5*(Y(I-1.J)*Y(I*1,J)) P(I,J)=.5*(P(I-1.J)*P( 1*1, J)) H8l./FL0AT(M-1) X(2,3)bO.

-93-

paQE

2,3) =-PES/M 3,?) = 0. 3,?) =-Y(2,3) 2.7) = 0. 2,7) = 0. 3,3) = .3333*(X(4,3)*X(3.2)*X(3.4) Y( 3,3) =.95*(Y(?,3)*Y(4.3)*Y(3.2)*Y{3.4)) P( 2,9) = 0. p( 2,3) =.5*(2-EPS) p( 3,3) =.5*(P(?.3)*P(4,3) p( 3.9) = 0. IP OLE =2*IP0LE-6 II = IPO LE-4 I3, II DO 33 CA LL F LIPZ(-1. Ill) S3 CA LL F LIPSK-I. Ill) sN*l 12 = M*2 II 0IP OLE DO 789 I- 2, 110.2 IN TFCs IPO LE-I-3 DO 789 J- 3, INTFC.9 X( I, J) S.5 {X( I, J-1 )*X{ I, J*l) Y( I, J) = .5 (Y( I. J-1 ) YC I. J*l) 789 CO NTIN UE II ObIP OLE DO 987 I- 3. 110.2 IN TFCs IPO LE-I-3 DO 987 J 2. INTFC y{ I.J) = .9 (X( I-l. J>*X( I*l, J) V( I.J) = .5 (Y( I-l. J)*Y( I+l, J) 987 CO NTIN UE RE TURN EN D
y( Y( X( X( y(
)
I

SUBROUTINE INTRR
DATA COMMON X(405,59).Y(405,52).P(405,52).M. CWY(40, 40), 72(40,40), CN, IPOLE.EPS. U. I?. ITER,ERX. IX, JX,ERY. lY, Jy.PRST, IS CITPSI,SR.GABY4,SRP.PABY4.C0N. I.J. ITPSTO.SBNR INTEGER EPS COMPLEX W.Z MleM+l SERYnO. ERXsO. SL0C1

JST

RES. H,

-94-

PAGE

lO

N2"N*N 0s.5*(Y(2,M)*Y(3.Ml) eRYABS(D-Y(2.Ml)) y(2,Ml)aD X(2.M1)0.

I2

INTrC5lP0LE-2
DO 111 JuINTrc.M 00 TO 1000

CONTINUE ini 111 CONTINUE

L0C2
I3,N INTrCsMAX0(2. IPOLF-I) DO 1 JelNTFCfMl 00 TO 1000 11 CONTINUE 1 CONTINUE IF (ITPSI) 3,2.3 21.22.21 2 IF (EPS)
DO 1 22 ERSI =0. DO 221 1=3.

!NTFC=MAX0(3. IPOLE-I)
DO 221 JbINTFCM 00 TO 2000

2211 CONTINUE 221 CONTINUE

00 TO 3 21 ERSIaO. DO 211 I3.N

INTFCeMAX0{3. IPOLF-I)
DO 211 JbINTFC.M WHYiiY( l.J)

A11./(Y( I*l,J)fWHY) A2sl./(Y(I-l,J)fWHY) A3al./(Y( I. J*1)*WHY) A4al./(Y( I. J-1>*WHY) A5ePABY4/(Al*A?*A3*A4) DsSRP*P( I,J)*A5*U1*P( 1*1. J)*A2*P( I-l. J)*43*P( EsABS(D>P( I.J) on TO 3001 IF (E.LT.ERSI ERSIiE
) )

T ,

J*l )*A4*P(

T ,

J-1

ISIal

JSIaJ 3001 P( I. J)sD 211 CONTINUE 3 CALL FLIPZ(-1. CALL FLIPZ(-1.

IPOLE-2) IPOLE-3)

-95-

PAGE

11

14

LOCaS MARKa IPOLE-3 DO 4 1=2. MARK INTFC =2*MARK-I DO 4 J=2. INTFC *-J-4) 4,4.1000 IF CONTI UE CONTI UE PSI) 5,6.5 IF LiPZd, IPni E-2) CALL LiPZd. IPnLE-3) CALL ITPSI ITPSI-1 RETUR LIP SI(-1. IPOLE-2) CALL LIP SI(-1. IPOLE-3) CALL
(

SMARKTP0L:

ITPSI
IF
(S2

ITP SIO

(F

DO 6?

61,69,61 S) la 3. MARK

INTFC 3*M ARK-I Ja 3. INTFC DO 62 QO TO 200 6?11 CONTI UE A21 CONTI UE 00 TO 7

61

MARK
DO 61

IP OLE-? In 2, MARK

INTFC xIP OLE-I Ja 2. INTFC DO 61 610,611,610 IF J-4 AlO WaCMP X( I-2)*H. J-2)*H) RES/W ZsCMP x(y I, J) Y( 1, J) Z2( I. )3 1./CaBS(7) = -AIMAQ(7)*Z2( I.J)**? WY( I. All CONTI UE MARK IPO LE-3 la 3. MARK DO 61 INTFC alP OLE-I Ja 3, INTFC DO 61 A12 P( I, J = 7 2( I, J)*P( 1, J) MARK IPO LE-4 la 3, MARK DO 61 INTFC MAR K-I*3 Ja 3. INTFC DO 61 WHYsW (I. J) Al = l. (WY 1*1. J)*WHY) A2 = l. (WY ( I-l, J)*WHY) A3B1. (WY I, J*1)*WHY) A4 = l. (WY I, J-1 )*WHY)
( 1 ) (
(

( (

-96-

PAGE

12

6121 613

614
7

inOO

A5=PABY4/(A1*A9*AS*A4) DaSRP*P( I, J)*A5*(A1*P( +1 J +A2*P (I -1 J)*A3*P( T,J*1)*A4*P( I. J-l) EABS(D-P( I> J) IF <e.LT ERSI) Gn TO 6121 SlSlBl SJSIaJ ERSI-E P( I,J)=D CONTINUE MARKbIPOLE-3 DO 614 I3.MARK INTFCalPOLE-I DO 614 J3. INTFC P( I, J)=P( I,J)/72( I.J) CALL FLIPZd* IPOLF-2) CALL FLiPZd. IPOLF-3) CALL FLIPSKI, IPni E-2) CALL FLIPSKl. IPnLE-3) RETURN DsSR*X(I,J)* GABY4*(X(I*1, J) X(I-l.Ji XM.J^I) X(T.,J-1)) EABS{D-X( 1/ J) IF (g.LT.ERX) SO TO 1001 ERXE
I
.

IXbI JX = J

inoi X( I, J)=D DsSR*Y( I, J)--GABY4*(Y( 1*1. J) +Y(I-1 . J)*V{ EsABS(D-Y( I, J) IF (E.LT.ERY) QO TO 1002 JYsJ ERYbE
lYal

J*1)*Y<

I.

J-l>

1002 Y( I,J)D 00 TO (ini.ll.l4>. LOC 2000 DsSRP*P( I, J)*PABY4*(P( I+l, J)*F(I-1 J)*P( EABS(D-P( I, J) GO TO 2001 IF (E.LT.ERSI ISlBl ERSIbE JSlBj 2001 P( I, J)5D QO TO (2001,2211.6211). LOC END
)

J*1)*P(

I.

J-H

-97-

APPENDIX II.

DIMENSIONAL PERTURBATION

In this appendix we report a correction to a

calculation described in reference [3],

Garabedian describes therein an entirely different

method for estimating certain parameters in axially symmetric

cavitational flows and jets.
In particular, the vena

contracta in three dimensions is treated as a perturbation

of the two-dimensional model with m as the perturbing

parameter, where m+2 denotes number of spacial dimensions,

as usual.

Hence the technique is know as dimensional

perturbation
The results lead to the conjecture that the ratio
of the radius of the jet at infinity to the radius of the

aperture, Y^^yY^,

may be expressed in m+2 dimensions as a

.

power series expanded around m =

Y(m)
_ Y(0) -

Thus

Y^^

Y^

2 _ - a^m + a2m +

Regarding m as a parameter in the equations for the vena

contracta, Garabedian is able to calculate
Yco(-l)

and, of course.

Furthermore, the derivative 8Y_/9m, evaluated

at
in

for constant Y^ = 1, is expressed as follows:

1

\
, [

^^0

= -

2
"^

3+a^
(1+a
)

^ tan'-'-a,. -l + ] tan ^

aya
^-!

db da

'

'

47ra +7Tb"+Yb

where
1/2
Y = 2 +

(1-a

- [4a +(l-a

'

+(l-a 2a+[4a 2,
o
,

;i

,,

2.2 -7rb/a, 1/2 e ) ]

^

a log

(l+a)"^

This integral
is quoted as

was evaluated numerically in

[3]

and the value

3m

= -

.650544

Interpolating the above data with a cubic polynomial in

5

= m/(m+2)

leads to the expression

^
^0

.6110 + .4857

.1110 6^ + .0143 6^

as an approximation to the power series.

6

For three dimensions,

is 1/3 and the result is

Y^TTT-

-"^^^^

leading to the estimate of the contraction coefficient

Y(l)
,

-99-

An indication of the accuracy of this estimate

is provided by the decrease in the magnitudes of the terms

of the polynomial.

The terms are, for m = 1,

.4857
=

.1619

.1110 6^ =

.0123

.0143 8^ =

.0005

and Garabedian arrives at an error estimate of 1/10 per cent.

Using our larger CDC 6600 computer we re-evaluated

the integral for 3Y-/3m
is
g

and found that a more accurate value

Y
,

= - .7072 + .0003 x-^ dm

leading to the cubic polynomial

.6110 + .5279

.1110 6^ - .0279 6^
0,

which fits the data at m =

1,

and

as before.

The revised computation of C

Yoo(l)

for m = 1 is then

C^ =
c

2
]

v l Yq(1)
/ \

.7737^ = .5985 + .0002

This, of course, is larger than the value in

[3]

and slightly closer to our evaluation, although it still

differs from the latter by .007.

A glance at the magnitudes

of the terns in the polynomial now reveals that the error

estimate m.ust also be revised upwards, for

-100-

.5279

.1760

.1110 6^ =

.0123
.0010
.

.0279 6^ =

Consideration of this data in the light of other

computations of C
leads us to conclude that the finite

difference, calculations are far Hiore reliable.

-101-

MAR 30

1971

DATE DUE

NYO-

C.2

Snider

N.Y.U. Courant Institute of

Mathematical Sciences
251 Mercer
St.

New

York, N. Y.

10012

SRe