The effect of air compressibility in the impact of a flat body upon a free surface

/

by

r.

THE EFFECT OF AIR COMPRESSIBILITY IN THE IMPACT OF A FLAT BODY

UPON A FREE SURFACE

Robert S. Johnson

Supported by

The Office of Naval Research

Report No. NA-66-8 Contract Nonr-3656(17)

September 1966

COLLEGE OF ENGINEERING

UNIVERSITY OF CALIFORNIA, Berkeley

THE EFFECT OF. AIR COMPRESSIBILITY IN THE IMPACT OF A FLAT BODY

UPON A FREE SURFACE

by Robert S. Johnson College of Engineering University of California Berkeley, California September 1966

It has been suggested that the high pressures exerted on the bottom of a ship's hull during slamming are developed in the air trapped between the hull and the water's surface. To

test this hypothesis, the two-dimensional, unsteady problem of the flow of air, where compressibility is accounted for, between a rigid, flat-bottomed block falling towards a rigid plane, is solved using a numerical method. The computed pressures exceeded those found experimentally by Maclean (Lewison and Maclean, 1966), and it is concluded that the deformation of the water's free surface must be accounted for in order to obtain agreement with the experiment.

To the author's knowledge, the numerical method, a modi-fied version of Sauer's method of Near Characteristics, is applied here for the first time and a maximum allowable time step, for this problem, is found by digital computer

CHAPTER II 11 CONTENTS Page ABSTRACT i CONTENTS ii NOMENCLATURE iv CHAPTER I INTRODUCTION i

Early Research into Slamming i

General Physical Model 2

Wedge Theories 4

Other Theories and Experiments 8

Specific Physical Môdel 14

STEP-WISE FORMULATION OF THE SLAMMING

PROBLEM 18 General Approach 18 Problem

I

19 Problem II 19 Problem III

21

Problem IV 23 Summary

24

CHAPTER III DESCRIPTION OF THE PROBLEM AND METHOD

OF SOLUTION 27

Problem Definition 27

Selection of Method of Solution

32

The Method

of

Near Characteristics 39

Numerical Method

44

Existence, Uniqueness, Stability,

and Convergence 53

Selection of the Mesh 56

Numerical Treatment of Boundary 59

Conditions The Corner 65 Shock Calculations 67 CHAPTER IV RESULTS 7:9 Generai 79 Principal Data 80

Time Step Limitation 82

Numerical. Results 83

Comparison with Experimental. Data 95

Effect of the Assumption of a

Rigid Free Surface 98.

CHAPTER V .SU1fl4ARY 101

Conclusions 101.

Suggestions for Future Research 102

Acknowledgements 104

APPENDIX.A IÑCORESSIBLE, INVISCID FLOW OF A. FLUID 105

BETWEE1N. AN INFINITELY WIDE,. RIGID BASE

AND A FALLING. RIGID. RECTANGULAR BLOCK

APPENDIX B TOTAL STRENGTH OF SOURCES ALONG B'C' IN 119

THE -PLANE (SEE APPENDIX A)

APPENDIXC DISCUSSION OF MAIN PROGRAMS (Including, .121

FORTRAN IV Listings)

APPENDIX D AUXILIARYPROGRAMS (FORTRAN IV Listings) 147

.APPEÑDIX E APPROXIMATE CALCULATION: OF FREE SURFACE 156

DEFORMATION

iv

NOMENCLATURE

a-

local speed of sound

speed of sound in undisturbed medium

a2, o,T

see equation (111.27)

/4 see equation (111.50)

b half-width of falling body

B

see equation (111.50)

C3

C4t/1 _{Jacobian elliptic function} Jacobian elliptic function

E

complete elliptic integral of the second kind acceleration due to gravity

ii height of body above bottom

I) J

indices for rows and columns in meh, repectively see equation (111.31)

see equation (A 7)

i

complete elliptic integrai of the first kind see equation (A 13)

f-

pressure

11- radial distance S entrophy

-4'7tfr1. Jacobian elliptic function

t time

respectively,

E

see equation (111.27)

V

velocity vector

V(t)

velocity of falling body

L4,r weight per unit area of the model

X. 7. horizontal and vertical coordinates,

' respectively, with origin at intersection

of body's centerline and rigid bottom mesh spacings in x- and y-directions,

respectively

location of shock wave formation

+ c_.

14 )

see Figure 10

complex coordinate Jacobi's zeta function

ratio of specific heats, taken equal to 1.4

JI

parameter

e density of fluid

qconstants

constant or time increment velocity potential

velocity of shock wave

coordinates in the

f

-plane complex coordinate

CHAPTER. I

INTRODUCTION

Early' Research into Slamming

In recent years there has been a steady trend toward schedule-keeping, as well as increased, ships' speed., even in heavy weather. This frequently, leads to viô1ent motions of the ships, which in. turn, lead to high oadings and consequent structural damage. The most violent of these motions lead

to bow emergence and. subsequent re-entry into the water, it

is generally believed that the impact phenomenon. at re-entry produces the high loadings and has, therefore, been the súb-ject of much. research, not only for ships (Schade et al.,

1964), but also for the similar problems involving seaplane floats .(Szebehely, 1952;. Chu and Abramson, 1961), planing boats '(J. D.. Pierson, 1950) and missiles (Moran, 1965). The

above references all have extensive bibliographies on the cited fields of research.

Although experimental. (Bottoniley, 1919) and theoretical (von' Karman, 1929) studies of the impact of seaplane floats date back to the 1920's, research in ship slamming did not generally appear until around 1950. There' ware earlier dis-cussions of the damage apparently aused by slamming by

Adam (1928), King (1934) and Kent.. (1948), but it was

Szebehely (1952) who first began to apply rational mechanics in the form of Wagner's (1932) theory to the problem.

General Physical Model

Before going into the various theoretical approaches, a schematic picture of the ship slamming phenomenon will be considered. Assuming that the problem is, in essence, only

two-dimensional, Figure 1 illustrates in an exaggerated way

the hydro-ae..ro-elastic problem. The entire structural field

*

of plating is assumed to be deñrmed by the built up air pressure. Likewise, the plating is assumed* to have additional deformation with respect to the over-all de-formation of the entire structural field. The water sur-face is also shown to respond to these air pressures.

Whether a significant amount of these responses takes place is still a moot point but this work intends to shed some light on these questions. In fact, whether the air's presence is significant needs substantiation and will be discussed below. Recent and on-going experimental work by Maclean (1966), Chuang (1965), Lewison and Maclean (1966)

and Glasfeld (1966) should provide much insight into these crucial questions. In this report a step-wise approach to the problem is set up which, in principle, should be able to take account of all of these possible complexities.

In constructing Figure 1, the phenomenon was considered to be two-dimensional. This could be implied from results of damage surveys reported by Schade et al. (1964). They state that "the extent of damage to the forward bottom structure

*

I

-1

-SIDE PLATING

BOTTOM PLATING SHIP LONG L FRA M) N G AIR WATER

FIGURE I.

A GENERAL

HYDRO-AERO-ELASTIC MODEL OF THE SLAMMING

PROBLEM

FOR A FALLING SHI P.

bottom surface; support structure is only affected in about l57 of the cases reported from some ship types and much less frequently from other types." If slamming damage occurs when the flat-bottomed part of the hull enters the water almost parallel to the water surface, then the location of the damage would tend to indicate a generally two-dimensional loading behavior. However, at the present "state-of-the-art", two-dimensionality is more hypothesis or conjecture than fact deduced from observations.

Wedge Theories

The earliest theoretical work was done by von Karman (1929) and Wagner (1932). They considered rigid, two-dimensional, wedge-shaped bodies impacting the surface of an incompressible, irrotational fluid at constant velocity

(see Figure 2). The presence of the air was ignored. Wagner's paper is somewhat compact, and J. D. Pierson

(1950) has worked it out in detail, extended it to large deadrise angles and applied it. Pierson indicates three possible approaches. In the first, the added mass and velocity field at the undisturbed position of the free surface are obtained by assuming that the wedge can be

FIGURE 2

WAVE SPRAY ROOT 2.c

-5-.

RIGID .WEDGE FREE SURFACE WATER

ALL WAVE. THEORY..

. CONTINUITY AT THE

FREE SURFACE.. GUARANTEES UNIQUE

SOLUTIONS.

...

RIGID WEDGE

FREE SURFACE WATER

FIGURE 2b.

_{ALL. SPRAY THEORY. CONTINUITY AT THE}

FREE SURFACE GUARANTEES UNIQUE

SOLUTIONS.

FREE SURFACE

WATER

FIGURE. 2c. WAVE AND SPRAY THEORY. CONTINUITY GIVES

INFINITE NUMBER OF SOLUTIONS, BUT

SIMILARITY AND AN EMPIRICAL SPRAY

THICKNESS RESULT IN UNIQUE

replaced by a body with the same beam as the wetted beam of the wedge moving in an infinite fluid. This results in a wave rising next to the body (see Figure 2a) and meeting the body at an oblique angle. The second is an assumption that all the displaced fluid goes into spray, and results in a picture like Figure 2b. The corner in the fluid ät the free surface in Figure 2b is pointed out by Pierson to be physically unreasonable. The third approach is a combination of the first two where an empirical spray thickness, g

at the nominal wetted half width, C , and a statement that

the flow be similar for each instant of time were required to obtain unique results. With respect to the first approach, the question as to what body to use in the infinite fluid to generate the up-wash field and added mass for the wedge was

left open above, Wagner and Pierson used a flat plate, whose width was treated as variable, and equal to the nominal

wetted width (see Figure 3a). Fabula (1957) has carried out the first approach, using all wave and no spray, for an

ellipse (Figure 3e). Bisplinghoff and Doherty (1952) have performed a similar solution for a diamond shape (Figure 3b). Meyerhoff (1965a) points out an error in the latter results

and corrects them. Recently Ferdinande (1966), apparently unaware of Pierson (1950), Fabula (1957), Bisplinghoff and Doherty (1952) and Meyerhof f (1965a), has repeated the solution for the diamond. However, his paper gives some details unavailable elsewhere. In a survey paper, Chu and Abramson (1961) have compared all of these wedge

-7-WEDGE

/FLAT PLATE

FIGURE 3a.

EQUIVALENT (EXPANDING) FLAT PLATE

USED TO OBTAIN UPWASH FIELD AND

ADDED MASS FOR WEDGE;

FIGURE 3b.

EQUIVALENT (EXPANDING) DIAMOND

USED TO OBTAIN UPWASH FIELD

AND ADDED MASS FOR WEDGE.

FIGURE 3c.

_{EQUIVALENT (EXPANDING) ELLIPSE}

USED TO OBTAIN UPWASH FIELD

theories. Three drawbacks are common to all the wedge theories:

They are not valid as the rise of keel goes to zero, i.e., as the wedge becomes flat bottomed

the solution breaks down. All the wave theories are valid for ttsmallhl deadrise angles although

Pierson has extended the combined wave and spray theory to large deadrise angles,

They predict very high values of the pressure, of the order of thirty times the stagnation pressure, in the neighborhood of the spray root for, say

a 200 deadrise wedge, and

. Fabula (1957) shows a graph which indicates that

these types of theories do not give even over-all total force results that agree well with expertnient for the early stages of the impact, i.e., for small

time.

It should be noted here that the flat bottoni case, that is, when the wedge angle goes. to zero, is probably of the

greatest interest for ship structures. Furthermore, the pressures at the early stages are of greatest interest and they are not predicted even approximately by these theories.

Other Theories and Experiments

Because of the drawbacks of the wedge theories, Egorov (1956) and Ogilvie (1963) formulated linearized versions of

the problem of the rigid body impacting on a compressible

fluid. These theories were carried out only for the flat

bottomed case (see Figure 4), but this is the case of greatest

interest. They result in a maximum pressure on the body,

equal to the acoustic pressure, where the accoustic pressure is defined as

Vo

where is the density of the fluid,

V is the velocity of the body, and

ct0 is the velocity of sound in the fluid.

Recent experiments by Chuang (1965) and Maclean (1966) have yielded pressure measurementS for flat-bottomed bodies, which were designed to be "structurally rigid" with respect

to the loads that acted on them, impacting on the water's surface. The tests were two-dimensional and were thereby related to Egorov's and Ogilvie's theories. The measured pressures were considerably lower than those predicted by

theory. In Ogilvie's and Egorov's theories the presence of

the air was ignored, i.e. , the body could he considered as falling through a vacuum. Thus, the experiments offer strong evidence that the presence of the air cannot be ignored.

Chuang (1965) also took high-speed motion pictures of his experiment and his photographs indicate what is apparently air trapped in the top layers of water immediately below the body just after the body makes contact with the water.

/77/

FIGURE 4.

"PHYSICAL"

_{PICTURE OF OGILVIE'S}

AND EGOROV'S

_MATHEMATICAL

-li-One would also expect that if the air's presence was ignorable then a flat-bottomed body would experience

pressures greater than or equal to those of a wedge-shaped

body. But Chuang's and Maclean's experiments indicate

the contrary, i.e., experimental pressures on flat-bottomed models are lower than theoretically predicted pressures on wedges. This is not to imply that experi-mental pressures on wedges would be necessarily less than experimental pressures on. flat-bottomed models.

Ex-perimental pressures on wedges might be higher than for f lat-bottomed bodies if the air has a cushioning effect, because the air could escape more easily from a V-shaped model thus cushioning it less than a flat-bottomed model.

Finally, a statement due to Fabula (1957) shôuld be

noted. In commenting on test results of two-dimensional

wedges of deadrise angles 100, 200, 30°,40°, and 500 reported by Bisplinghoff and Doherty (1952), Fabula notes that theoretically predicted and experimental

deceleration time histories are in good agreement except for the 10° deadrise model. At 10° the theory predicts decelerations that are too high. Furthermore, Bisplinghoff and Doherty (1952) give sketches (Figure 8 of their paper)

of the water surface, determined from high-speed motion pictures, which show that for the 10° wedge the water

surface is depressed from the undisturbed location just away from the body! For the other wedge angles the water surface was always above its undisturbed location. Fabula mentions then in a footnote "Since their record for

= loo shows an unusual depression of the free surface, air pressures created by the flat wedge may have

in-fluenced the case appreciably. ti

All the above theories and experiments were for structurally rigid bodies, Although they imply that the presence of the air seems quite important, the flexibility of the structure may be equally important. Fortunately, a giant step forward in this area has been taken by

Meyerhof f (1965a,b). He has extended Wagner's (1932) flat-plate theory to the case of a wedge-shaped body with

flexible bottom structure. As indicated by Diagram 27 in Meyerhoff (l965b), the pressures predicted at early times for a flexible-bottomed wedge-shaped body are quite close to and slightly less than those predicted using Wagner's (1932) theory. At later time the results for

a flexible bottom oscillate about those for a rigid one. Thus it seems that the introduction of structural

flexibility into the prediction of the pressures can only be a partial explanation, at most, because Wagner's results are poor for early time.

200

o

-SYMBOL EXPERIMENTER MODEL (PSF)

DESCRIPTION FLAT BOTTOM

U-. SHAPE

FLAT BOTTOM FLAT BOTTOM FLAT BOTTOM

u. Li

Ii

A

V

s

A

V

s

A

V

.

A

U

V

A

V

-U CHUANG 70 U CHUANG 210

-V

MACLEAN 274

w

A

MACLEAN

354

U MACLEAN

484

0.2

0.4

0.6 0.8 LO 2.0 4.0 6.0 DRQP HEIGHT (ft)

FIGURE 5. EXPERI'MEN1LLY DETERMINED.

for various flat-bottomed bodies. The U-shaped section reported by Chuang (1965) had a flat bottom of one-third the area of the flat-bottomed model tested by Chuang

(1965). If one assumes that both models were about the

same weight, this implies that the U-shaped model had about three times the weight per unit area of the f lat-bottomed model. This tends to confirm the data, due to Maclean ILewison and Maclean (1966)], that higher weight per unit area results in higher pressures. Some of Glasfeld's recent work (Glasfeld, 1966) also supports this conclusion.

Lewison in Lewison and Maclean (1966) gives a "quasi-two-dimensional" theory for the problem of a rigid body

falling through the air towards the water's surface. He calculates the pressures exerted by the air on the body

and is able to get remarkable agreement with Maclean's experimental data. His maximum pressures occur before the body makes contact with the water. Unfortunately, there are some somewhat arbitrary parameters in Lewison's theory which must be determined from experiments or from more exact theories.

Specific Physical Model

Chuang (1966) suggests that the maximum pressures do not occur until the body contacts the water but that some

-15-of the air is t.rapped between the body and the water and is forced into the water0 This he contends reduces both the speed of sound and the densi.ty of the fluid so that the acoustic pressure is now considerably lower. In

light of all this evidence it seems very reasonable to conclude that the presence of the air must be included

i.n a successful theory. But both Lewison!s and Chuang!s

points of view are defensible and need a definitive experiment to determine which, if either, is correct. Figure 6 illustrates the situation. The numbered points represent times at which the body might make contact with the water. If the contact is made at points 3 or

4 then an approach in the spirit of Lewisons is

appropriate. If contact is made ät point 1. then

some-thing like Chuangs hypothesis is indicated. And if contact is made at a point near 2, then the problem becomes horribly complicated, involving a matching of two solutions. For the purposes of this paper it will be assumed .that contact is made at a point, near 3 or 4, hopefully nearer 4.

Thus it wil,l be assumed in the following that the problem is one of gas. dynamics. The body and the water surface will, be viewed as boundary conditions, The

general approach will be to solve the exact. equations of motion in two dimensions using numerical methods,

loo

AT Ç. OF FLAT BOTTOMED

_{BODY. DROP HEIGHT:5f1.}

0.55

TIME (sec)

FIGURE 6.

TYPICAL

EXPERIMENTAL PRESSURE

HISTORY

AT A SINGLE PRESSURE GUAGE.

-17-The "state-of-the.-ar.ttt af..numeric.al methods

.for..high-speed digital. computers. can tréat either two-dimensional, unsteady flow of viscous, incompressible fluids. or

tWo-dimensional, ..unsteady flow of inviscid., compressible

fluids. Since this. is. obviously a. compressible fluid

problem, it will be necessary to assume that the. fluid

is inviscid...This is not to say. that viscous effects cati be ignored. in. setting up. the. mathematical model. 'A

well-known example of. such .an effect is .the' Kutta-.

Joukowski condition for the trailing, edge' of wings in inviscid fluids. Further comments on this subject wi'l.l be' made, as appropriate, in the main, text.

CHAPTER. II

STEP-WISE FORMULATION OF THE SLAMMING PROBLEM

General Approach

This Chapter will show how to formulate the slamming problem in a series of steps that will, hopefully, lead to better and better approximations to physical reality. The step-wise approach is based on the idea that different pieces of the same mathematical model have different response times. If it is true that the response of one part of the system is completed before another part can respond at all, then these two parts can be effectively decoupled in the calculations. A second idea that will be used is that it should be possible, where response of a part is "small", to neglect the response

of the part when computing forces acting on it. Then, the response of that part of the system can be solved for sep-arately under the now given forces. This response, in turn, can be fed back into the original calculation to find a new set of forces acting on the piece that, hopefully, will be close to those originally calculated there. If the system were a very simple one, this second idea would not need to be invoked, because all the pieces of the system, and their interac:tion,could be handled simultaneously. As will be

seen later, even the most simplified model to be handled herein approaches the capacity of a large-scale, high-speed

-19-Problem I

Figure 7 is the most simplified model of the problem

illustrated by Figure. 1.0 In Figure 7, the model and the

water's surface. are taken as rigid, and the. air is assumed

inviscid. The rigidity assumptions amount to a hypothesis

that the response times of t.he bottom plating, and the water's surface., are riuch greater than t.he duration during which the main pressure hump acts (see Figure 6) The

neglection of the air's viscosity is somewhat arbitrary and, as mentioned i.n Chapter I, due to the limitations of the current. "state-of-the-art." in high-speed computation0 The

corners of the model were given special treatment in the actual computation, and this will, be discussed in more detail in Chapter IL Whether this is a reasonable model

is open to question0 Current experiments are expected to provide empirical data which will confirm or deny the

rigidit.y assumptions. But., whet.her they are confirmed or

denied, Problem I., as outlined in Figure 7, i.s the first step

in the problem solving sequence.

Problem II

Figure 8 illustrates the second step, Problem II. Here,

it is assumed more likely that the water will. respond during the slam than that the bottom plating.wii.10 The deformation of the water's surface 'is obtained from a given function.

RIGID SURFACE

FIGURE 7.

PROBLEM I.

PERFECT GAS

(AI R)

-21-The function is derived separately by solving the linearized potential-theory problem for the water with the pressure distribution obtained from Problem I given at the free

sur-face. The. solution to this potential-theory problem is

given by Wehausen and Laitone (1960) and can be used in either of two ways

L

The water-surface deformations and ve.locities are given a priori for all time, or

2. . A formula for the deformations and velocities

as a function of time and,, say, average pressure below the body at the location of the undisturbed water surface is given, .and as new average

pressures are computed during Problem II, the formula i.s applied to find the water!s motion. This phase, Le.., Problem II and the potential-theory solution, can be repeated iteratively until predicted pressures and pressures used to find the wat.er!s motion'

agree to within some acceptable tolerance.

Problem III

Figure 9 illustrates the thir4 step, Problem III. This is analogous t.o Problem II in' that the bottom plating's res-ponse is computed separately. Pressures from Problem II would be used to predict the form of deformation, and then some function of average pressure and time would be used

WATER

PROBABLY SOME

SIMPLE FUNCTION,

e.g. SIMPIL COSINE OR PARABOLA.

AIR

-i----APPROXIMATION BASED ON

PROBLEM I PRESSURE

DISTRI-BUTIONS AND POTENTIAL THEORY

FIGURE 8. PROBLEM .

RIGID

BO DY

AIR

APPROXIMATION BASEDON RESULTS

BLEMS I AND

FIGURE 9.

_{PROBLEM m.}

-23-in actually calculat-23-ing the bottomts response dur-23-ing Problem

III. It is anticipated that deformations of the plating

will stay within 1tsmall-deflection" plate theory,. although using "large-deflectio&' plate theory is not out of the question. For simplicity, probably something like simple supports at the ends and no intermediate supports would be

used. Some of the experimental work alluded to above may

suggest whether Problem III is even necessary, although for completeness, it might still be carried out. The structure above the bottom is regarded as rigid because in actual ships the slamming loads occur locally and the overall structure responds much later, if at all.

Problem IV

Finally, Problem IV is to calculate the response of the bottom structure for all time. It is expected that even if the bottom plating experiences deflections during the slam, it will experience considerably larger defor-mations after the slamming loads are over. That is,

accel-erations due to the slamming loads will cause these defor-mations. For stiff enough bottom structure, the.pressure oscillations. (see Figure 6) subsequent to the main hump have little effect on the plating's response. This last statement is somewhat contrary to Meyerhoff's (1965 b) results, where plating oscillations .caused alternate reinforcement and relaxation of pressures as against the

rigid body case. However, this current work is a completely different approach from Meyerhoff's, and the verification of the statement will only come about through an experiment where a flat, flexible bottom is tested. Thus the loads for the bottom plating can be obtained from Problem III (or II, or I, if that is all that is required) and then a completely separate and uncoupled plate-response problem can be carried out, using those loads as input.

Summary

What is being assumed here is that a very complicated

hydro-aero-elastic problem can be reduced to one in gas dynamics, where only the significant, and anticipatedly small, effects

due to the elasticity of the bottom and deformation of the waterts surface, are included. It is anticipated that the approach will predict the main hump in Figure 6, up to at

least point 3. The rest of the hump is not so critical,

and it could be supplied by extrapolation based on comparison with experimental data.. Having obtained the hump, no more pressure calculations need be done, and the prob.lem is reduced to obtaining the time-dependent structural response to the given time-dependent loading.

The structural problem may not be trivial, but it is certainly considerably easier than the originally posed hydro-aero-eLastic problem illustrated by Figure 1. Since ships experience permanent deformation of their bottom

-25-plating, we expect to consider plastic effects. Fortunately, there exists considerable guidance for this part of the analysis. Keil (1961) presented a very, simple total-energy approach

that has met with considerable success for plates subjected to blast loads. It assumes that the plate acts as a membrane and no plastic yièld occurs. until the entire field of the plate goes plastic instantaneously.. Wah (1960, 1961) and Nagai (1962. a, b, e, 1963) have treated. the infinitely long

plate usinga.plastic-hinge type of analysis. 'f necessary, a numerical approach could be used where the plate is

divided into finite elements. In the numerical approach rauher arbitrary stress-strain curves could by incorporated into the algorithm. However, current. finite-element techniques could not handle an elastic perfectly plastic matérial, but could handle a strain-hardening mteria1.

While the above approach is no.t guaranteed to lead to. a completely successful solution to the slamming problem, it does have the merit that it breaks the problem into discrete steps, all. of which seem solvable within the bounds of present knowledge, . Of course, if any of the hypotheses about. response

times are .bad.l.. in error, the. approach may fail. And the

important implicit assumption, that. the.peak.pressures. on the body occur before. the body makes. contact with the water, may also prove to be incorrect.. The truth of. .these assumptions

should be partiálly validated by the results of the current calculations. Final confirmation rests with experimental

evidence.

Problem I is the key to all the succeeding steps. If

a numerical solution to it can be found, then Problems II and III, if they are needed, will simply be small variations of Problem I. The balance of, the present paper will be

devoted to the solution of Problem I. Problem I is the

necessary first step, as well as being of intrinsic interest on its own merits.

Momentum:

2L2L

*2r2L.

-7L2,

-27-CHAPTER III

DESCRIPTION OF THE PROBLEM AND METHOD OF SOLUTION

Problem Definition

The problem to be.. solved is that of an infinitely high, two-dimensional, rigid block falling towards an infinitely wide, rigid plane (see Figure 7). The block is assumed to fall, under t.he influence of gravity but the pressure build up beneath the block will be taken into accòunt in determining the velocity history of the block. The fluid will be assumed inviscid and the f low will be considered isentropic, The

governing equations* .are, therefore, the following (Stanyuko-vitch, l96O)

Continuity:

(IIL2a)

(IIL2b)

*Subscripts

denote partial differentiation with respect to the subscripted variable,

where is given by initial conditions on

_p

and Assuming that air is a perfect gas and introducing the

definitiòn of the speed of sound, O..

then by suitable manipulations, one finds

where _è' = 1.4.

Further manipulations of equations (111.1) through (111.5) yield the equations of motion in terms of

U, 2/

and Q..

only, [see, e.g., Holt (1963 b)1:

o,

(III.6a)

(III.6b) State:

-29-I

(zc +-,)=

o,

(III..6c)

This form of the equations of motion was selected so that in the numerical calculations only Q. and not böth

and would have to be carried along. If only one of

or

r

were carried along, then equation (111.3) would

constantly have to be used to recover the other as needed. Equations (11L6) are hyperbolic (Courant & Friedrichs, 1948) quasi-linear partial differential equations. Quasi-linear means that they are non-Quasi-linear but that the highest derivatives of the dependent variables occur only to the first power, although the coefficients of the highest deriva-tive may contain powers of lower dervaderiva-tives. The problem itself is a mixed boundary-value initial-value problem. If

one introduces the speed of sound in the undisturbed fluid, the initial conditions at every point in the flow field are

O)

_(11107a,b)

The boundary conditions are (see Figure 7): Along Ä:

u-

= 0, = arbitrary, O.. = arbitrary. Along : 2h. = arbitrary,

1T

velocity of the falling body*, = arbitrary. Along

b:

24.. = 0,

if

= arbitrary, O-. = arbitrary, 0, = 0, 0. Along : 24. = arbitrary, ?.= 0, = arbitrary,

14=

0, 0, 0. (III. 10) (111.11)

*The _{body's velocity is taken as -p with a correction for} the pressure build-up between the body and rigid bottom.

-31-Furthermore, equations (111.7) must be satisfied along some curve connecting a point, on to a point on [see Figure

71.

This curve represents the front of the disturbed flow

field. Observe that the boundary conditions on 24. and

2P

in equations (111.8), (111.9) and (111.11) are just slip conditions at solid boundaries. The conditions on various derivatives come from the symmetry and antisyminetry of the various parts of the flow field. They may be more easily visualized by referring to the equivalent problem pictured in Figure A2 in Appendix A,

The use of "slip conditions" on and amounts to an assumption that the respective boundary layers on

and are very thin compared to the height of the body above the bottom, When the body gets near the bottom this may not be true, especially if the flow is laminar. By

making a cru4e incompressible calculation it can be seen that the Reynoid's number of the flow is above the critical

Reynold's number for a pipe of the same diameter as the distance from the bottom to the body, when the latter distance is about one inch. Thi.s would imply a turbulent

flow during. t.he time period of greate.st interest. However,

it has been pointed out to the author that this entire phenomenon is over so quickly that a turbulent flow may not have time to develop. This is an open question and might prove to be an interesting research problem.

Observe also that along

b ?J

O. This means that on

the centerline the flow is not the one-dimensional piston

problem. U. O reflects the fact that there is flow

away from the centerline, just off the centerline, and that this flow relieves the pressure at the centerline.

One of the first steps in approaching this problem was to search the literature for the corresponding in-compressible-fluid problem. Such a problem had not been worked out. It was accordingly solved and since it was not trivial has been included herein as Appendix A. The complex velocity is given for an arbitrary point in the flow field as equation (A. 17). Unfortunately, a numerical quadrature

remains to be carried out, but since the solution was of no direct use, it was left hanging at that stage.

Selection of Method of Solution

The question now arises of how to solve the problem represented by equations (111.6) through (111.11). The

non-linearity of the governing equations puts the problem out of reach of an exact analytical solution. A

linear-ization is not appropriate because the problem is of interest precisely when the dependent variables have undergone large

changes. Since there are characteristic dimensions,

simi-larity methods are also inapplicable. Thus numerical methods were turned to.

-33-The numerical solution of partial differential equations received a great impetus during World War II from such people as Southwell (1946) who developed relaxation methods to such

a high degree.. Such techniques were limited to the amount

of computation a man was willing to perform. With the advent of the modern high-speed digital. computer this restriction was dropped in favor of hcw fast. the machine one was working with was and how large a computing budget the investigator

had. During the 1950's both the mathematics and computers evolved to the extent that two-dimensional (one space and

one time or two space) problems, even for non-linear equations, were routine. Not. only were answers obtained but existence,

uniqueness, stability and convergence of the methods had been proved (Forsyt.he and Wasow, 1960). These calculations

were carried out on the first generation (see Lapidus, 1962 for the definitions of computer generations) of high-speed digital computers of which the IBM 650 and IBM 704 were

typical. About 1960 a second genèration of computers,

such as the Remington Rand LARC and IBM 7090., became operational. Concurrently, interest in doing

three-dimensional problems developed (Forsythe and Wasow, 1960; Fox., 1962)... Although nearing its end, the second generation of computers is still with us. Much computation and theo-retical work has been done on three-dimensional problems. Groups in Australia (Richardson in Alder et 1964),

England (Butler, 1960), Germany (Bruhn and Haack, 1958),

Russia (Godunov, Zabrodin and Prokopov, 1961; Belotserkovskii and Chushkin in Holt, 1965) and the United States (Alder, Fernback and Rotenberg, 1964; Lax and Wendrof f, 1964; many others) have worked on these problems during this period. Unfortunately, mathematical definitions and results governing stability and convergence have been limited to linear equations or linésrized versions of the non-linear equations.

There are at least four basic approaches to obtain

difference equations that are approximations to the differential equations (111.6):

Finite-Difference methods - The derivatives in equations (111.6) are directly replaced by finite-difference approximations to the derivatives.

Characteristicmethods - Real characteristics

exist for hyperbolic partial differential equations. Equations (111.6) can be put in an equivalent form in which derivatives are taken along space curves called characteristics. These derivatives are then replaced by finite differences.

Conservation methods - The equations of motion

are put into conservative form and then differenced. For example, Noh (in Alder et al., 1964) takes

momentum and density as dependent variables. Equations (111.12) are in conservative form and

-.35-(III3) when the indicated vector products are carried out.

where V

= (i

?i).

Equations (1.iI12) all involve the operation

where

f

is a scalar0

The Conservation methods then develop and use

difference approximations to the operation indicated by (11L13)0

4. Integral-Relations methods - These methods (developed

in Russia) have just recently been published in English by Beiotserkovski,i and Chushkin (in Holt,

1965).. The equations Of motion are written in divergence fórm and then integrated along strips in one direction.

(ru)j

(ev)

(III. 12)

±

(,u,,

LÇ

,

= Fi

J (111.14)

where q

,

/

, are given functions.

Equations (111.14) are in divergence form for two independent variables ) and . The integration

in one direction leads to a system of ordinary differential equations which then can be handled by the well developed methods of the numerical solution of ordinary differential equations.

Furthermore, these methods can be formulated in Eulerian or Lagrangian coordinate systems or in combined

Eulerian-Lagrangian coordinate systems.

The choice of method is difficult. Certain generali-zations can be employed to rule out certain classes of

methods. For instance, the. Integral-Relations method was un-known to this author when the current problem was being

pro-grarnmed. Lagrangian methods initially label a set of

particles and then follow their history. Since most of, if not all of, the initially labelled particles of some regular grid beneath the body would have passed out from under the body by the time the time period of greatest

sparsely defined. In this problem the flow field beneath the body is of primary interest, and thus Lagrangian coordinate systems were ruled out. However, Problems II and III would probably be simplified by the use of combined

Eulerian-Lagrangian coordinate systems.. Lagrangian coordinates could be used near the water's free surface and the body's bottom and an Eul.erian grid t.o define the rest. of the flow field. Beyond these statements the choice is somewhat arbitrary.

Before proceeding further, the basic idea used in solving two-space dimensional., unsteady problems should be mentioned. One is given all values of the dependent variables (and their necessary derivatives) at some time, t0. Then one has an equation or equations of the form of equation (111.15).

(111.15)

To advance _4. from the given time, to, to a new time,

to +t, one approximates

Eby

using its arguments at.

time t0, replaces

24

by the approximation

and arrives at

(t0+t)='zL(z')+ ¿Z

(III. 17)

The derivatives appearing as arguments of are obtained

approximately by numerical differentiation. Equation (111.17) could be improved iteratively by replacing the arguments of

F

by average* values. This is called an implicit method, where the values of the function at time

t0 +t play a

role in determining themselves. When only values of the

dependent variables at time t0 are used to calculate values

at t0 +t, then the method is called explicit. An explicit method is usually the first step in an implicit method.

The method of Characteristics was chosen for the com-putations. It was believed that for the same mesh a

Characteristics method would be more accurate than a first-order Finite-Difference method. The Lax-Wendrof f second order Finite Difference method is stable in some norms and unstable in others for linear systems. By example, Burstein (1964, 1965) has shown some situations where it is unstable and so it, too, was ruled out. The most promising. Conservation method looked to be that due to Noh (in Alder et al., 1964). A conversation with Noh indicated that his method took two programmer-years to program and should be avoided because of this.

*For _{example, replace 2Iby}

where

2r(1)

t+t)

is the latest estimate for that quantity.

-39-The Method of Near Characteristics

Holt (l963a, 1963b) has surveyed the literature with respect to Characteristics methods and has concluded that Sauer's method of Near Characteristics (Sauer, 1963) is the

most promising. Furt.hermore, to this author's knowledge, nobody has ever tried to carry out a calculation using Sauer's method. Therefore., Sauer's method was chosen for the calculations. For completeness a discussion and descrip-tion of Sauer's method will be included.

A disadvantage of early three-dimensional Characteristics methods [e0g0, Butler (1960)1 was that a complicated

inter-polating scheme was required to find the location of points lying back along a characteristic. Sauer's method simplified this by f in.ding three independent characteristics lying in one plane in the three-dimensional. space0 To cover the spáce a grid can be lai.d out on the (x,y)-plane and the

charac-teristics can be found in (x,t)-planes Figure 10 illustrates this situation, In the equations to be obtained, partial

derivatives in the y-direction will appear0 These are approximated by numerical differentiation using values in the parallel (x,t)-planes.

To find the governing equations in characteristic form, equations (1.1.1,06) are multiplied, respectively, by the

three constants and The resulting equations are summed and the terms involving y-derivatives are moved

u(J,I+I) > O

u(J,I)<O

u(J,I-I) <u(J,I)

- -

CHARACTERISTIC CURVE

J-I

J

IJ+I

FIGURE IO.

GEOMETRIC

_{INTERPRETATION OF}

-41.

to the right-hand side.

-cÇ2r

-,-/

where

G3

=

a

Or, rearranging the left-hand side,

Lc* +(a,u )j

[cjì

+

We now ask if there are

cr, C. ,

0

_{such that the}

derivatives i.n square brackets are all in the same direction, Or, are the following ratios equal.

aa

r

p

cC3t2.

where i.s a constant?

Equations (I.IIl.9) can be rewritten in matrix form as

or,

2L±c2,

T3=ZL.

u--r

O

c3Ct

o

_o

Q

C

₃

These equations possess a solution if and only if the deter-minant of the coefficient matrix equals zero. Hence,

(u-2-)[(2L- 2)

o

I'

7:

>(1

Oj

ci

_3/j

-ri

>(1, o,---)

1,

0)3

0Z]

_o,

=- 0.

(111.20) (111.21)

There will be a set of ( O , , G1 ). corresponding to

each . They can be obtained by putting equations (111.21)

one at a time into equations (111.19).

Putting (111.22) into equation (111.18) results in three independent equations:

-)

= _Yzrzro2,

)=iì,*o

(lIli 23)

__

= - 2J__

-

ct

where

=

-(f)EÌ*

2.

Equations (III 23a) and (III23b) are valid along curves governed by the differential equations

(111.24)

cLx

and equation (111.23e) is valid along

2A,

(III. 25c)

Equations (11L23) are three equations in the three unknowns

ZL

,

ZI

and &. They are implemented numerically by means

of a scheme similar to that described by equations (111.15) through (111.17).

Numerical Method

Having developed the above method, Sauer (1963) indicated that the next step was to replace the derivatives appearirg in the differential equations by finite-difference

approxima-tions. Unfortunately, as will be seen later, he did not try

the method out.

The first step in developing the numerical approximation. to equations (111.23) is to locate the curves along which they are valid. These are characteristic curves passing through point (J, I, t0 +t) in Figure lO. The second step is to find the intersection of the characteristics with the

t0-plane. This is done by approximating equations (111.25). Equation (III.25a), governing the "+" characteristic, will be used far illustrative purposes. The differential

equation is replaced by a difference approximation

-45-.

4'

,.

where

U

and a. are average values along the "+" characteristic.

Then i? and are replaced by

U

and . at ( ')', t0),

respectively. In an explicit method that is ali that is done In an implicit method, once an estimate for

U

and Q. is

calculated at

t0 +t, then thevalues of Uand

_Q. at times .4

t0 and t0 +t are averaged to replace 'Z2and Q. . Since

everything i.s known in équation (111.26) but

Xt,

one can solve for it, If one does the samé thing for equations

(III.25b) and (III25c), all the intersections of the characteristics with the t0-plane are found. Let us define

z)-

Q

These areall found by linear interpolation between the appropriate mesh points. This is the chief feature of Sauerts method, Interpolation was only required in the x-direction,

Now equations (111.23) are differenced, Equation (III.23a) will be used as an illustration:

where % is

_(J-i)A

is

(r-I)4

and replaced by

¿L

(

z

/

z)

- u

(J) Z/)

)

d (111.29)

Similar approximations are used for

and Q,

Ex-pression (11.I.29) is a central-difference approximation to the derivatives and makes an error of

_0(4/).

Doing the same thing to equations (III23b) and (III.23c) and solving for 2L(J, I, t0 +4t), 2 and

c2,

one obtains

u(

=u-)-(-zr

o-

(j,

r

-t)=

-+2,_a 7

_¿J,

(111.30)

7r(J

1, t+ht,)=

?1

/,

caJ.

M(I.t0t)

ztt)

=

W'+

_(111.28)

-47-Again, all the terms on the right-hand sides are known at time, to, and so values of the dependent variables can be

computed at time, t0 +t. In an implicit method more

accurate approximations to all the variables on the right-hand sides can be introduced after the first explicit cal-culation0 However, it was found that the values of

and were essentiall.y unaffected by iteration and that values of derivatives did not change significantly over the time period t0 Thus an explicit method was. used in the calculations reported herein. Normally difference methods of this sort, subject to time-step requirements governed by stability considerations, would be well behaved away from the boundaries of the flow field and the only problem re-maining would be the correct handling of the boundary con-ditions0 This was not the case for the equations (111.30). The procedure described above for advancing one time step does not work0

The failure of the algorithm manifests itsel.f by spatial oscillations in the computer results. Much effort was

invested in determining that these oscillations were not due to instability due to taking. too. large a time step. The key to discovering this was that the oscillations were of a period of about (see Figure lIa) in the y-direction but smooth in the x-direction.

u,v

or a

O CALCULATED POINTS

POSSIBLE SOLUTION THROUGH POIN S.

FIGURE lia.

TYPICAL CALCULATED RESULTS

USING SAUER'S METHOD WITHOUT

AVERAGING IN y DIRECTION.

TRUE SOLUTION

-- AVERAGE SOLUTION

s

FIGURE lIb.

THE TRUE SOLUTION CANNOT BE

REPRESENTED BY A GRID OF

SPACING

_{S, BUT ONLY SOME}

-49-understood, but such a variatiOn in the computed results can be attributed to the inevitable errors due to the difference approximations to the derivatives and to the rounding error in the computer. However, their persistence can be explained. The central-differencé approximation to derivatives, formula

(111.29), by which gradients in the y-direction are taken into account, just i.s not sensitive to any gradient of about

I period

2.4f

or less. The major contribution to this

insensitivity is that the value of the function at the point at which the derivat.ive is being taken is not involved in calculating the derivative. In the x-direction the value of the function at the point being calculated plays a major role in determining thé values of the function back along the characteristic curves. This serves to average out any random

f luctuatiáns The explanation of why the error in using

formula (11129) is not is that the derivation of this error estimate presupposed that the third and higher derivatives of the function were of O(l) In.Figure lla a possible solution corresponding to the computed function values was drawn in and this shows that the third and higher

dérivatives need not be small. That is, the use of formula (11129) presupposes that the function is always smooth in the y-direction0 A theoretical insight results from all of

this If the true solution to a problem is as pictured in Figure llb, then the grid shown in Figure llb could not

could represent would be some average solution. Therefore, choice of a grid immediately ±u.les out finding any components of the true solution higher than some frequency which is re-lated to the mesh width. Hence, one can only claim that a numerical solution represents the average behavior of the

true solution. Arguments that the physical situation requires smooth behavior are alsó suspect because one needs only think of a Krmn vortex street behind a moving cylinder to find a highly oscillatory, asymmetric, real flow where one would normally assume some sort of smooth, nionotonic, symmetric behavior. The only conclusive check for high frequency

com-ponents would be further grid refinement, but usually this is too expeisive in computer time.

Because of the above, the choice of frequency components one will let into the computed results is somewhat at the investigatorts discretion. He must certainly screen out oscillations such as those illustrated in Figure lia. Since all the high frequency terms are already screened out by

the mesh, nothing is essentially lost by screening out a few more.

The author developed a means of accomplishing this screening out without any loss in accuracy of the method. The procedure desc.ribed.above had errors at. least as large

as LJVj/J.

What one does LS t.o use .an approximation to the

%I*

-

O

H

-51-example, one calculates z.( I + 1,

t0)

and

'LL. (

I - i,

t)

when one is interpolating for

[E24.(

t,

I,

t)].

Expanding as Taylor series,

I)

t0).i-y z--'(r)t)

-I- 7

t0),

_E

(i

A,4/

ZL (

i /

r)

=

z

(x. '

I

-

I) t,)

"X)?&) z),

E(Er-L74(r-fl4/).

Adding A211to both sides and dividing by , one gets

¿(I1) t0)

_ut+2«I-L

₎₌

O(í

(III3l)

A choice of 2 amounts to averaging, by means of the trapezoidal rule and

4.

= 4 corresponds to Simpsons first

rule. As can be seen by examining the program listing in

Appendix C, this complication to the programming adds about 307e to the calculation time at a typical mesh point. There is a theoretical interpretation to this averaging. Let us consider the region of dependency at time, to, of a point at (J, I, t.

+t)0

It is approxirnatel.y a circle of radius

with center at (J, I, t0). If not too much variation in the values of the dependent variables occurs in the x-direction, then Figure lia may be used to visualize the situation. The

average value of a dependent variable over the interval

[I -

½, i +

½] is less than the value at row I (the reverse

is true at row I + 1). A calculation scheme should somehow take this into account because the value of a dependent

variable at (J, I, t0 +4t) depends on values of the dependent variables in its entire region of dependency. The pure

central-difference approximation doesn't do this but the technique represented by the formula (111.31) does, In

general, was taken as 2, However, near some flow-field boundaries where y-derivatives were known to be zero, was

taken as 4 because this simulated the zero slope there. At least one other investigator has found similar problems. Trulio (in Alder et al., 1964). states that "Centered schemes of differencing can lead to severe non-physical oscillations," In light of the above discussion about Krmn vortex streets, the author questions the use of the word nonphysical. However, whether physical or non-physical, the severe oscillations were present. Trulio used a rather elaborate identification scheme to find such oscillations and then, having found them, used a technique to smear them out.

-53-suggested a least-squares application of hi.s method, In two space and one time dimensions, this amounts to carrying out his method twice with characteristic planes taken in the x- and y-directions and averaging the results. He suggested that this might improve the methods accuracy. i.t would also

probably eliminate the hi.gh frequency fluctuations because finding the values of the dependent variables back along characteristics in both directions would introduce an averaging into the computation in both directions,

Existence, Uniqueness, Stability and Convergence

This section will attempt to clarify the mathematical foundations upon which the computations are based. _Courant and Hilbert. (1965) point out that for first order, quasi-linear hyperbolic systems of equations t.he solutions exist and are unique0 This is, of course, for the ar1ytical problem and not the numerical approximation. A necessary condition for uniqueness of the numerical scheme is that the domain of dependence of the difference equations shoul.d include

the domain of dependence of the differential equations. This immediately leads to a limitation on the allowable time

step to

Lt

't

Sufficient conditions for the stability of the numerical method and the convergence of the numerical solution to the

analytic solution, as the mesh width is decreased, are still open questions for all numerical schemes that have been

applied te quasi-linear systems of equations. The most that has been done is to assume that the flow field is quasi-stationary over the time step

¿t

and, hence, that the coefficients of all derivatives are constant. This reduces the stability analysis from a non-linear problem to a linear problem. Unfortunately, the results of such an analysis are only valid for local time. Burstein (1965) has shown, by a computational example, that such an analysis of the Lax-Wendroff Finite Difference method was not sufficient and, in

fact, that the method was unstable under certain conditions. He attributed this to the neglect of the non-linearities.

Despite the drawbacks of the linearized approach, there is no other guidance. From a linearized approach to a Conser-vation method, Godunov, Zbrodin and Prokopov (1961) give an

allowable time step as

(111.33) where

f?v% (a-!/)

_rj

and

-Z)

-55-For arid t! and 2f small compared to O-, this becomes

,L4::jLC

o-Noh (in Adler et aL, 1964) states that if

e

everywhere in the flow field, then

4

f»1A,

(zX)Lf)

- (111.35)

He states, "For the general case the following restriction has been found satisfactory in practice:

LVtZ

-J ZA.

I +

/X-J1tI

(111.34)

(III, .36)

He was also dealing with a Conservation method, yet inequality (111.36) reduces to

L_1-

,1

for

and Uand1T very small compared to

O...

In the light of the foregoing and because it was not the purpose of this paper to go into the mathematical analysis of the method, it was decided to incorporate an arbitrary parameter, called CFLCON, in the computational scheme, where

¿t = C FLCON

rrr-T)J

(111.37)

CFLCON was set initially at 0.5 and was capable of external adjustment as the computation proceeded.

Selection of the Mesh

in almost any elaborate numerical calculation many implicit assumptions are built in. These range from when the algorithm should regard a computed number as zero to the assumption that a situation is at standard temperature and pressure because a particular value of the density of a

substance is part of the algorithm. Such implicit assump-tions are usually undesirable, but as any person experienced in programming knows, they are difficult to identify and keep out.

In this problem the choice of mesh carries with it many implicit assumptions. Many of the decisions about the mesh had to be made without any experience with the solution

because the mesh must be determined before a solution can be obtained. And once one has invested a considerable amount of money in computer time one is usually reluctant to repeat all calculations because a slightly different mesh looks as though it might be better. That is, the solution affects the choice of mesh and the choice of mesh affects the solution. One hopes that this interaction will not be excessive.

-57-Very soon after beginning this problem it became apparent that the problem's running time would be on the order of from 4 to 10 hours on the available iBM 7040/7094

DCS. This meant that an efficient algorithm with the least

number of mesh pouts for the desired resolution was

manda-tory. Immediate]", this led to a grid with variable spacings

because this would allow detail where detail was required and general trcnJ where they were judged sufficient. The

single greatest obstacle to the successful algorithm was the treatment of the corner To the author's knowledge there is. no discussion of corners in subsonic compressible flow any-where in .the analytical or numerical literature. Many

numerical problems have been solved about cornered shapes but the authors never discuss this problem The treatment of the corner largely determined the character of the grid. A full discussion of how the corner was handled will be

included in a succeeding section of this chapter.

Figure 12 depicts the grid at 0.12 seconds into the problem. Besides the fixed grid there is a movable row, labeled MB, which is tangent to the body's bottom.. Although it moves, dependent variables calculated at it are found by using equations. in Eu.lerian coordinates. Necessary informa-tion at old time is supplied by interpolainforma-tion before each time step. The fixed grid itself is not permanent. The

body is allowed to fall from row 1.1 to r.ow l0. When it reaches row 10 the grid is changed in the y-direction.. The overall height of grid i.s kept fixed; the intervals below

COLUMNS 9

un"..

I1HIII1UI

_A

°° !!IIIHIIIUU

i.uuuuuui'i'..

CORNER DETAILS 5. 062 5 FIGURE 12.

THE MESH BEFORE

THE FIRST REGRIDDING IN THE Y- DIRECTION. COLUMNS:

/

13/21 / 312/ 36////4í///// // ///46 / / //// / /// //// / /

// /

ROWS: 21 17 SI COLUMNS SPACING If?) I-3 0.25 3-II 0.50 Il-21 0.25 21-31 0.50 31-36 1.00 36-41 2.00 41-46 4.00 46-51 8.00 IO II 13 15 I? IS 21 22 23 f ROWS 13 12 _{MB IO} 9 8 13 MB 6

-59-row 13 are reduced to 0.9 their previous value, and the intervals above row 13 are stretched to equally subdivide the distance from row 13 to row 21. This results in a new grid where row MB starts out again at row 11. A good

resolution is thereby maintained in the region of primary interest below the body.

Nümerical Treatment of Boundar Conditions

As far as possible the boundary conditions were applied to the algorithm in a natural manner. Along the boundaries

and the analytical bouñdary conditions were

directly applied in the algorithm, e.g., on fwas taken as the body's velocity and on

29.

was set equal to zero. Along , the algorithm was reformulated taking the

charac-teristics in the y-direction. [Incidentally, permutation of variables doesn't do this easily because signs change. The

safest course is rederiving the characteristic equations

from equations (1116)] This caused

_Ç

and O

to appear

in the equations and then they were set equal to zero

Along the highest row (row 21) and the column furthest from the centerline (column 51) special treatment was

neces-sary. The region of dependency of the points on these

boundaries covers part of the mesh but also extends beyond the mesh out to the edge of the undisturbed flow field. This is schematically illustrated by Figure 13 which, unlike the algorithm used, shows y-characteristics However, no values

of dependent variables are available to the right of row 21. This presents two possibilities:

Extend the mesh to encompass the entire disturbed flow field, or

Develop some way of transferring the boundary conditions at the edge of the undisturbed flow field back to the boundaries of the mesh0

The first possibility is very unattractive because the body falls to the bottom in about 0.56 seconds and this means the grid would have to be about 600 feet square x 0.56

600 ft.]. Even with a coarse mesh beyond the grId shown in Figure 12, the resulting number of mesh points would be

pro-hibitively large. The second possibility is very attractive, but the procedure finally adopted was somewhat arbitrary.

The acoustic approximation to a cylindrical piston expanding at constant velocity was carried out. That is,

the surface velocity of the piston was regarded as small compared to the speed of sound and the equations of motion were linearized using their ratio as a small parameter0 The resulting equations possessed a similarity solution. This

solution related the speeds of sound and radial velocities of two points on the same radial line

I ₂

Ii

._I'tli

V,

ct0t

-61-COLUMN J

u CHARACTERISTIC

\

"-" CHARACTERISTIC

,,, I'<

s,

-

_-r

REGION OF DEPENDENCY

OF(J,2I,t)

/

UNDI STURBED FLOW 4 t

-Lv-4 7 18 9 20 21 ROWS

FIGURE 13.

REGION OF DEPENDENCY OF

DIFFERENCE EQUATIONS AT

POINT (J,21

,

t).

and

ct0t

i ---"

()21

(f"')

1cttJ

a2

(111.39)

where ,'L is the radius.

Radial lines were drawn from the initial location of the corner and the formulas (111.38) and (111.39) were used to extrapolate values at row 20 and column 50 to row 21 and column 51, respectively. While this is an admittedly approximate procedure it serves the practical purpose of relating the variables at the outer boundary of the dis-turbed flow field [at a distance o_0t from the initial position of the body's corner] to the variables at the mesh boundary.

The variable position of row MB and the change in vertical mesh width at row 13 caused some problems in numerical differentiation. Since the central difference approximation (11L29) is accurate to O(y2)it is natural to look for a scheme of numerical differentiation that will maintain that accuracy. If one fits a parabola through the

three points 7J(J, I + 1), ?J(J, I), and

2f(J,

I - 1)

and then differentiates the parabola at point

(J,

I), the result is an approximation of accuracy The

situation is even worse than this would seem to imply, as can be seen by reference to Figure 14. The "true" function

-63-SS

/-APPROX 1MAI I NG PARABOLA ROWS

FIGURE 14.

EX4MPLE OF BAD APPROXIMATION

_TO

"TRUE" SOLUTION

_{THROUGH VU ,12),}

V(J,13) AND V(JI4).

14

"TRUE"

SOLUT ION

4J

¡5

through computed points 1, 2 and 3. Then the "true" function decays away gradually over computed points 3, 4 and 5. If one tries to approximate at row 13 by the derivative of the parabola passing through points 2, 3 and 4, the result can be catastrophically poor. If, however, the function is behaving very smoothly [almost linearlyl any reasonable scheme will give a good approximation to the derivative0 Thus at rows MB and 10 below the body, where everything was very smooth, parabolas were fitted and differentiated to find approxi-mations to the derivative of the function. More details of this will be found in Appendix C0 The solution to the problem at row 13, where large gradients could exist, was to

inter-polate for the slope of the function rather than approximating the function and differentiating the approximation0 Since the vertical mesh spacings are constant above and below row 13, the derivatives can be found accurate to at rows 12,

14 and 15. The

Au's

are different, but that only introduces

a constant which can be absorbed into the

Q(

₎ notation. Unfortunately, the theoretical error is still but dealing with the slopes of the functions puts one on an in-herently smoother curve, in this instance, and the practical consequence is a generally smooth variation of the function across row 13. Hence, practically speaking, better approxi-mations result,

-65-The Corner

Next to the oscillations discussed in the previous section the treatment of the corne.r was the, most difficult aspect of developing a successful algorithm. Initial attempts to approach the problem directl.y and calculate what was going on at the corner met with failure0 All these attempts re-gar.ded the corner either as square or of infinitesimal radius0 The gradients in the dependent variables, particularly the pressure, were extreme, despite the fact that the fluid was not treated as incompressible, and the finite mesh could not cope with them so that large destabilizing oscillations

would result.

Af ter studying the grids used in handling obstacles with square corners, where the fluid was assumed incompressible or compressible but moving subsonically, (Noh in Alder et aL,

1964; Harlow and Welch, 1966) it was observed that these grids avoided the corner. That is, the location of the córner was taken so that no mesh point landed on the corner and so the corner's direct treatment could be avoided0 Because of these previous successes, it was decided to emulate this procedure

here0 In Figure 12 'it canj be seen that the corner was located

one-quarter of a mesh width to the right of column 13 The

justification behind avoidance of the corner is that one is

not really interest.ed in what happens there, as long as, the rest of the flow field is computed correctly0 _{Apparentlythe}