Numerical study of large amplitude free surface motions

THE PHYSICS OF FLUIDS

TCE U'

'

VOLUME 9, NUMBER S

2bOlOTiUOf9aS.

Scheepshydron!eth8flpa

Numerical Study of Large-Amplitude Free-Surface MotionP

ekeIwe 2, 2628 CD Deft

Fauwis H. Hai.ow AND J. EDDIE Wztcx O1 788? F8t O1518182

Los Alamos Scientific Laboratory, University of California, Loa Alamos, New

Mexico-(Received 5August 1965)

-The marker and cell method for high-speed computer was used to investigate the motion of a fluid whose free surface has been perturbed by an impulsive sinusoidal pressure. For gravity pointing ont of the fluid, the resulting motion exhibits Rayleigh-Taylor instability, whose progress is followed

from low amplitude to the asymptotic bubble and spike stage. The effects of small and large viscosity

are contrasted. For gravity pointing into the fluid, the large-amplitude oscillatory motion exhibits a shift of resonance frequency and other effects of nonlinearity. Comparisons of results are presented wherever previous analytical or experimental work was available, and many new, aspects of these types of flow are discussed. The calculations are based upon the full, time-dependent, nonIin

Navier-Stokes equations for an incompressible fluid.

TRODUCTION

-WHEN

impulsive force on its free surface, the naturea box of fluid has been subjected to an

of the resulting motion depends upon whether the body acceleration (gravity) points downward (into

the fluid) or upward (out of the fluid). In the former casè, the motion is stable, while in the latter it exhibits the well-known Rayleigh-Taylor instability.

We have investigated both types of motion using the marker and cell technique' to solve with

high-speed computer the full. nonlinear, time-dependent

Navier-Stokes equations of motion for a viscous

incompressible fluid with a free surface. The examples

to be discussed are for motions in two space dimen-sions, in a" rectangular box. The method. is equally

applicable to three-dimensional studies, that are precluded at present only by the limitations of the

available computers. . .

There are several purposes which computer studies such as these fulfilL First, they serve to test directly. the applicability of the equations on which the

com-putations are based, by. direct comparisons of the.

'results with physical experiments. Second, they give information about the processesfor which laboratory

experiments are too difficult or costly to perform. Third, and perhaps most important, they resemble accurately controlled and recorded experiments, furnishing data upon which analytical models can be based and tested. Thus, computer studies com-plement the work of the laboratory experimenter and of the analyticál theoretician, but certainly do

not replace eithr. ' .

All calculations reported in this paper were

per-formed on the IBM '7030 (Stretch) Computer. The

'F. H. Barlow and J. E. Welch, Phys. fluids 8, 2182

(1965).

842

figures shown'g particle configurations, velocity

vec-tors, .or préssure contours were processed directly

from the computer by. the' Strómberg-Carlson

SC-4O) Microfilm Recorder and are not retouched or

otherwise altred.'

,

THE NUMERICAL METHOD

- The inarkei: Wand cell technique has been described

in detail in Ref. 1, so -that only .a brief discussion

.of the basic features will be presented here

The region through which the fluid moves is

divided into a fixed (Eulerian) coordinate mesh of

rectangular cells, each .of dimensions öz by ¿y. This enables the fluid to be represented by a finite1

set of values of the appropriate, field variables.

In the marker and cell method, these variables are primarily the pressure and the components of ve-locity. A typical calculation will utilize about two

thousand such cells.

- The fluid itself is represented by a set of marker

particles. The motions of these are determined by a weighted average of the velocity components of the cells through which they move. In addition to show-ing the changes of fluid, configuration, the particles

enable the surface position to be followed, so that

the computer can automatically incorporate the

free-surface boundary conditions. The particles do

not otherwise enter into the calculations.

In the solution of the equations of motion, the

cell dimensions give the finite-difference intervals for approximating the space derivatives. In

addi-tiön, the time derivatives are similarly approximated, so that the' calculation proceeds from the initial

conditions through a 'sequence of cycles, each of

finite duration Ot. The finite-difference approxima-tions are derived from the Navier-Stokes equaapproxima-tions written in the following conservative form:

au

i+i;+ +g,+vs+),

8u3 du ô fa2u ô2u\

8z' duv f't9V a2v\

--I-±j, +=g,±+r),

öu öv

ô+ôyO

The velocity components u and y are parallel,

respec-tively, to the z (horizontal) and y (vertical) axes. The ratio of pressure to (constant) density is

des-gnated by ç; y is the kinematic viscosity coefficient; and g,, g, are the components of body acceleration gravity). The only approximations involved in the

slútion of these equations are those coming from

the finite-difference representation. In particular, all of the nonlinear terms and viscous. effects are fully

incorporated.

The calculational steps through one cycle of time

advancement are as follows: (1) At the beginning of the cycle, all necessary cell-wise data and

par-ticle coòrdinates are available, either as initial

con-ditions or as rempining from the previous cycle

The velocities satisfy the finite-difference analogy

to the incompressibility equation. (2) Pressures are calculated The Poisson's equation for deter-mining the pressures is based upon the requirement

that the new velocities calculated with them will

also satisfy the incompressibility equation. The finite-difference Poisson's equation is solved by an iterative procedure. (3) The new velocities are calculated from the finite-difference versions of the momentum equations. (4). The marker particles are

moved to new positions. Their velocities are

cal-culated by a linear interpolation among the four

nearest cell velocities.

These four steps complete the basic parts of a cycle of time advancement. Actually, there are

-many minor aspects, presented previously,' which

have been omitted from this brief description.

Questions of accuracy of a computing method

such as this one can never be satisfactorily answered except by abundant comparison of the results with

those of experiments or other solution techniques.

Certainly it is necessary that the mesh of

com-putational cells be sufficiently fine to resolve the essential details of the flow field. Also, the time

step per cycle (öt) must be cim1l enough to ensure that field-variable changes per cycle are almost

everywhere much less than the values of the

vari-ables. (We find, in fact, that if the configuration

of the marker particles is plotted every cycle and the results assembled into a moving picture, then

IS sufficiently small if the projected viewing

ap-pears smooth. Ordinarily,

this means that any

marker particle must move somewhat less than a

cell width in. a cycle.)

The time increment per cycle is limited also by

considerations of calculátional stability.

Hydro-dynamic stability is often of interest to investigate; the full Navier-Stokes equations are completely

adequate for the study of many of the types likely

to occur in nature. Written in finite-difference form,

however, the equations contain the additional pos-sibility of numerical instabilities, of which the in;.

vestigator must always be wary. For the marker

and cell calculations, this consideratioñ means a

further restriction upon the value, of öt, whose

stringency is increased with increasing viscosity. The matter is considered in more detail in Ref. 1.

The viscosity enters in another way into the computational aspects of any marker and cell method

solution. Some Eulerian computational techniques

require nonzero viscosity in order to maintain stability We have found, however, that marker

and cell calculations can be performed with

con-siderable accuracy whenz' = O. me results in most

respects are indistinguishable from those obtained by using a small, but nonzero viscosity coefficient. The only observed difference is in the smoothness of the

marker particle configurations. For z' = O, the

ir-regular results may appear "turbulént," while 'with

small viscosity there is a much "cleaner" appearance to the particle arrangement. We have, therefore, generally incorporated viscosity into "nonviscous" studies and considered that z' is sufficiently small if

the boundary-layer thickness is everywhere much smaller than the cell size. In such cases a free-slip

boundary condition is then used at every rigid walL

Appreciable effects of true viscosity can also be

incorporated into the calculations; one extreme case

is discussed in this paper. .In such examples, the

marker and cell technique results aré particularly

smooth, provided that öt is sufficiently small. Rigid

walls then require, of course, a no-slip boundary

condition if the expected boundary-layer thickness becomes at least comparable to cell size.

-As previously discussed,1.2 the marker and cell method has been shown to possess a wide range of

applicability to fluid dynamics calculations. Thé

results shown in this paper give only a small sample. RAYLEIGH-TAYLOR 1ISTABILITY

With gravity pointing upward, out of the fluid,' the free-surface irregularities are unstable. Thus,

'F. H. Harlow, J. P. Shannon, and J. E. Welch, Science

F. H.. HARLOW AND J. E. WELCH

after an initial sinusoidal pressure impulse oñ the

surface, the fluid continues to fall upward with iii-creasing amplitude of the surface wave. Gradually

the shape of the perturbation itself changes as

higher harmonics couple in through the nonlinear effects. The leading parts narrow into spikes while

the trailing parts broaden into bubbles.

-This describes qualitatively a special case of the

Rayleigh-Taylor instability, a fluid dynamics effect

'which has received particular attention since the

war, by Taylor,3 Lewis,4 and more recently by

Birkhoff5'°. and Chandrasekhar,7 who summarize much of the earlier work. These and other authors

have described the grOwth of the instability through

the low-amplitude (linear) phase, for which

ana-lytical techniques are readily available. Such

mod-ifying. features as surface tension, viscosity, ¿nd

continuos layering have also been considered.

Experiments show, however, that when the

per-turbatioù amplitude exceeds about four-tenths of the wavelength, the linear-theory predictions are no

longer correct. After. a complicated transitional stage, a spike-and-bubble pattern is reached, in which, ideally, .the spikes fall freely and the bubbles move with constant speed. Analytical studies of

this final phase have proved, somewhat successful in their. predictions. of spike-and-bubble motions, but the transitional phase has been analyzed.only

with simplified models' (see Fermi's study, as

pre-sented by Birkhoff'), by series expansions of limited validity Mnge (see Emmons, Chang, and Watson8), or by numerical techniques. Some of the numerical mèthods : for incompressible fluids have been sum-marized !by Birkhoff,'° but none have been ex-tensive 4ough for useful analysis.

The p1esent technique,.while applicable to only

onè fluic with a free surface; enables computer

calculatidns to study in detail the full development

of the flow for both low and high viscosity. Accurate results have been obtained with relatively little

expenditure of computer time, so that we are able

to show the effects of several parameter variations.

The units of all calculations have been scaléd in

such a way that the magnitude of the body

ac-celeration is unity and the half wavelength of the I. Taylor, Proc. Roy: Soc. A201, 192 (1950).. D. J. Lewis, Proc. Roy. Soc. A202, 81(1950).

5G. Birkhoff, Los Alamos Scientific Laboratory' Report

LA-1862 (1955).

'G. Birkhoff, Los Alamos Scientiflë Laboratory Report

LA-1927 (1956).

S. Chandrasekhar, Hydrodynamic and Hydromagne&

Stability (Oxford University Press, London, 1961).

W. Emmons, C. T. Chang, and B. C. Watson, J.

Fluid Mech. 7, 177 (1960). '

pressure pulse is L = 4.8 distance units (a number

which was convenient for irrelevant computer usage

reasons). Units of the ratio of pressure to density, denoted by p/p = , and of the time t, are thereby

uniquely established. In every example, the pressure pulse was applied to the plane surface of a quiescent fluid. This pressure pulse is given by

çc3=A6(coskz,

V

in which k = ir/4.8. In the above units, .A = 5.0 provides a high-amplitude impulse, causing the

fluid to "coast" quickly to the nonlinear phase of

motion, while A = 0.5 enabled additionally the analysis of the early linear phase. Computer studies' were made, for both values of A, as well as for two

values, of D, the initial depth, and for two values

of the kinematic viscosity. Finally, one calculation was performed to compare with the analysis of Emmons, Chang, and Watson.8

The results of one of the computer runs are shown.

in Figs. l-5. For this example, A = 5.0 (large

impulse), D = 8.0 (deep water), and i' = 0.01

(negligible viscous effect). Figure 1 shows the marker

particle positions at a sequence of times,

ifiustrat-ing the growth of instability well into the

spike-and-bubble phase. The sides of the confining box

are walls of reflective symmetry; for computer

efficiency,' only the minimum permissible region

was 'covered in the calculation. Figure 2 illustrates another type of output available from the computer.

Each line is a velocity vector whose origin is the

center of a computational cell, whose length is proportional to the lòcal fluid speed, and whose directiOn is the local flow direction.

-The lines in Fig. 3 are isobare, which show the

dramatic changes occurring in the pressures at early times. A peculiarity, of the plot program

slightly shifted the isobara relative to the framing

lines of the box. The obvious horizontal correction puts them into their proper orientation.

An analysis of the linearized equations shows that the interface displacement Z should follow the behavior

Z = A(/gD)'sinh [t(gD)1] coskx,

(1)

inwhich'.

,

V

kDtanh(kD). (2)

For early times, V

' VV' V

Z' = * (At/D) ces kx, . (3) referred to as the initial impulse motion. '

V

81)plles only to the earliest part of the motion. In fact, as shown in Fig. 4, the effects of nonlinearity

become pronounced well before there is significant

difference between the predictions of Eqs. (1) and (3).

The solid curves of Fig. 4 give the surface displace-ment histories of the spike and the bubble. The line splitting the two at early times shows the initiai

impulse prediction of Eq. (3). The dashed curves of Fig. 4 are the spike and bubble speeds. It may

be noted that after the time t = 0.4, the spike speed increases at a uniform rate. The result from the calculation is that the spike tip falls freely with exactly the constant gravitational acceleration, g =

+ 1.0, as expected.

The bubble velocity has been replotted in Fig. 5,

in several ways. Also shown is the radius of

cur-FIG. 1. Fluid configuration in Rayleigh-Taylor instability for the calculation with A _{= 5.0, D = 8.0. The times are} 0.00, 0.25, 0.50, 0.75, 1.00, 1.25.

F1o. 2. _{Velocity vectors corresponding to the frames}

shown in Fig. _{1. A vector whose horizonal component is the}

width of the box indicates horizonal component of speed 4.8?

with all other vectors in proportion. The vertical component scale has been cut in hail for_{ease of visualization.}

vature of the bubble. Birkhoff5'° has discussed

studies of several investigators concerning _bubble motion, and jointly with Carter6'9 has developed

an analytical prediction for steady-state bubble speed. The predicted result° is indicated _{by one of}

the horizontal lines in Fig. 5. A similar result by Garabedia&° is also indicated in the figure, and

our agreement with the latter is excellent.

Lewis,4 on the other hand, measured somewhat

greater bubble speeds, but Birkhoff5 suggests that

his accelerations were considerably _{greater than}

re-G. Birkhoff and D. Carter, J. Math. & Mccli. _{6, 769}

(1957).

'° P. R. _{Garabedian, Proc. Roy. Soc. (London) A241, 423}

846

F. H. HARLOW AND J. E. WELCH

Fia. 3. Isobar plots for the calculat on shown in Fig. 1, at the times 0.1, 0.2, 0.3, 0.4, 0.5, 0.6. Interval between isobars is 0.5, with the zero isobar at the surface not plotted. The internal isobar at ç = Ois indicated by a mark at the side.

ported, thereby making bubble-speed analysis

un-certain.

The radius of curvature of the bubble in our calculation, measured in units of perturbation

wave-length, is R/2L = 0.39, a value slightly higher than

that indicated by Birkhoff,5 who assumes the value

0.35 and who also mentions the unpublished

ex-perimental result of Duff, R/2L = 0.26. The Duff measurement may be a misprint; the data for his

results are not available. Reference to Lewis' paper, in which interface photographs _{are reproduced,} shows a variation in R/2L from 0.35 to 0.40, which

results should be independent of any uncertainty in the acceleration of his system. We are therefore confident in the correctness of the computer value.

To further verify the accuracy of our calculations, we next examined a case similar to the one de.. scribed above, but in which A = 0.5, so that the early, linear phase could be compared with Eq. (1). The calculation is of interest only in that it showed excellent agreement with the analytical predictions. The results described so far are for a fluid whose

depth is sufficient to be considered infinite. For a shallow fluid, the behavior is expected to be

some-what different, particularly for the bubble. Figure 6 confirms this expectation for a fluid with initial

depth D = 2.0, whose conditions were otherwise identical to those for the calculation described in

Figs. 1-5. A detailed comparison with Fig. 4 shows that both the bubble and spike motions are impeded. With bubble amplitude limited to an asymptotic value of 2.0, the constant-velocity phase can never occur. The spike, however, quickly reaches the free-fail stage, as expected.

For these calculations with negligible viscosity, the accuracy of the results has been assessed

prin-cipally by comparisons with theory and experiment8 at the earliest or latest phases of the motion. For

the intermediate transitional phase, there is very

little comparative material. One exception is

fur-nished by the calculations of Emmons, Chang, and

Watson,8 whose theoretical treatment included

ef-fects to third order in the amplitude. Graphs of

their interface configuration show clearly the be-ginning of the spike-and-bubble tendency, with a large lateral shift of the null-amplitude point. To compare with one of their results, we performed a calculation starting with a perturbed interface and

no initial impulse. The results are in close agreement with those of Emmons, Chang, and Watson, with a

small discrepancy arising at the latest times that

they show. By that time, it may be supposed that

-TIME

Fia. 4. Amplitude (solid) and velocity (dashed) of spike and bubble as functions of time for the calculation shown

0.8 06 II' 0.4 o o, z w 6 0.2 (5) for which an approximate value is ' = 0.91262. Equation (4), whose validity depends upon the approximation of linearity, shows a decelerative

coasting of the fluid asymptotically to a deformed

configuration. The initial impulse motion is the same as that of Eq. (3). A comparison between Eqs. (1) and (4) shows that at the earliest times after

the impulse, the viséous damping must dominate

over the gravitational acceleration.

More generally, the full eigenvalue equation for

the exponential growth rate coefficient can be written

1t2-4ar+4a2t[1+(r/a)1'l} =0,

(6)

in which and a are dimensionless measures of the

time coefficient and of the viscosity, respectively,

n2g'k1,

a2

kYg'.

(7)

[The amplitude is proportional to a weighted sum of terms in exp (ni).] Of the two roots, one is real

and positive, while the other is complex, with

relatively large negative real part. An analysis of

the equation shows that viscous effects are negligible if a2 « 1; Le., if 2

« gk3. Under that circumstance, the appropriate solutions of Eq. (6) are =

i

-2a + -2a1 añd =

1 - 2a - 2ia, with fairly

good accuracy for a 0.1. For large values of a, the real root approaches (2a1, while the complex

TIME

Fia. 5. Bubble radius and velocity as functions of time for

the calculation shown in Fig. 1, compared with the predictions of Ref. 6and IO.

their neglect of higher-order effects contributes to

an inaccuracy which, in comparison with the

com-puter results, shows a slightly narrower spike and

broader bubble.

It is well known that the effects of viscosity and surface tension on RayleighTaylor instability can

be considerable.7 In the initial linear growth phase, the effect is to severely limit the growth rate of short wavelength perturbations, with the result

that a dominant mode is determined. In the later

phases, the shapes of both the spikes and the bubbles

may be significantly altered..

As an extension of the idealized flows described above, we have calculated several in which viscosity

plays a strongly modifying role. (A technique for

including surface tension in the calculations has been proposed, but not yet tested.) With the surface

impulse initial condition, there is immediate com-petition between two energy sources. The gravita-tional body force accelerates the flow as before. The viscous dissipative forces, however, tend to

slow down the flow in a msrnner which can be de-scribed, for the extreme case in which g = 0, by the surface behavior

Z = (A cos kz/1a4')[1 - exp ('kt't)]. (4)

The dimensionl nimìber is the nonzero root of the equation

Fia. 6. Amplitudes (solid) of spike and bubble as functions of time for calculation with A = 5.0, D = 2.0. Datum points show scatter in calculated spike-height positions.

848

F. H. H.ARLOW AND J. E. WELCH

Fio. 7.. Cnflgurations and velocities for the calculation which otherwise is the same as Figs. i an 2. The scales are the same as in the times aire 1.0,2.0, 2.5.

-one is prbportionai to a, with negative real and imaginary parts...

The sequence of fluid configurations illustrated

in Fig. 7 demonstrates the strongly modifying effects

that large viscosity can have. The conditions are

identical to those of the problem ifiustrated in Figs. 1-5, except that the coefficient of viscosity is y = 2.0. Correspondingly, a = 1.058, and the real root of Eq. (6) .j = 0.377. [In contrast, the low-viscosity

calculations, with y = 0.01, have a = 0.0053 with

only a 1.1% effect on. the roots of E. (6).]

A comparison between Figs.. .1 and. 7 shows,

first of all, .the expected viscous effect of slowing the motion considerably. An elapsed time of 1.0 in Fig. 1 corresponds, in spike motion, to a time of 2.5 in

Fig.

_{7. In addition, there are other qualitative}

high-viscosity

illustrated in

those figures;.

differences shown by the comparison.. Viscosity in-creases the width of the spike and narrows the

bubble, so that the basic harmonic of the

perturba-tion persists longer. .

Comparisofi with Fig. 2 shows that there also

are qualitative differences in the appearance of the field of velocities. Strong motions of convergence nto the spike characterize the low-viscosity case, while the large-viscosity effect is to make the ve-1 cities near the spike center more nearly vertical.. A quantitative comparison of the viscous effect is visible in Figs. 4 and S. (Note the necessary dif-ference in time scale between the two.) In Fig. 4

the theoretical initial impulse solulion lies between

the calculational results for bubble and spike,

whereas Fig. 8 shows that the strong effect Of vis..

cosity is sufficient to retard both spike and bubble

motioñ to well below the initial-impulse curve. Also shown in Fig. 8 is a plot of Eq. (4); the close

agree-ment with the ihean of spike and bubble motions

indicates the small effect that the gravitational

body force produces. The effect, however, is not

negligible. At the latest times, the bubble speed

is given by u5(2Lg 0.184, while its radius of.

curvature is given by R/2L = 0.3 13. These, which

are nearly the final asymptotic values, may be com-pared with the low-viscosity results, 0.24 and 0.39,. respectively. .. .. .

As before, we wished to examine in much more detail the early, linear phases óf the process. Thus, this same large-viscosity run was repeated with a

much smdler initial impulse, A = 0.5. Figure 9 shows the full history of the configuration

develop-ment, while Fig. 10 shows quantitatively the spike-.

and-bubble motion. The calculation was carried far enough for the steady state to be established.

The final spike motion in Fig. 10 corresponds the acceleration g = 1.0, while the bubble is

char-u o I--J Q.

a

I 15 TIME

Fia. 8. Amplitudes of spike aùd bubble as functions of time for the calculation shown in Fig. 7, to be compared with

acterized by u8(2LgY1 = 0.181, and R/2L = 0.303. (Jlese values confirm the statement that the cal-cuintion shown in Fig. 8 had nearly reached its steady state.)

NONLINEAR OSCILLATIONS

Theoretical calculations1' and confirming ex-riments'2 have both indicated that the resonance frequency of a box of oscillating fluid changes with amplitude in a manner which depends upon the

mean fluid depth. For a deep fluid, the frequency

decreases with increasing amplitude, but if the depth

Fia. 9. Fluid anfiguration for the high-viscosity

calcu-lationwithA =O.5,D =8.O.Timesare 1.5, 3.0,4.5,6.0,_{7.5, 8.75.} 11L Tadjbakhsh and J. B. Keller, J. Fluid Mech. 8, 442

(1960).

D. Fultz, J. Fluid Mech. 13, 193 (1962).

4 6 8

TIME

Fia. 10. Amplitudes of spike and bubble as functions of time

for the calculation shown in Fig. 9.

to wavelength ratio is sufficiently small, there is an increase in frequency.

To investigate this effect, we changed the

pre-viously-described calculations only insofar as direct-ing the body acceleration into the fluid. All scaldirect-ing and dimensions remained unchanged; _{i.e., L = 4.8} (half wavelength), g _{= 1.0, y = 0.01 (low viscosity} in every case). The only variations _{came from}

changing the depth D and the

_{strength of the}

initial impulse A. Thus, each calculation _{started from}

the half-wave impulsive loading of the fluid surface

and followed the motion through half a period of oscillation. A typical result is shown _{in Fig. 11,} in which are plotted the crest and _{trough motions} as functions of time. In the example, D _{= 4.0 and}

A =1.0.

Fia. 11. Amplitudes of spike and bubble as functions of time fortheoscillatingflowinwhjchA = 1.0,D = 4.0.

F. H. HARLOW AND J. E. WELCH

The results of linear analysis predict a surface

motion

Z =

A(/gD) sin [t(g/D)J cos.kx, (8)

which also is plotted for comparison in Fig. il.

Just as in the Rayleigh-Taylor instability, one of the principal manifestations of nonlinearity is seen in the enhanced amplitude of the spike side, and

the reduced amplitude of the bubble side. (Through

the second half of the period, the spike and bubble

sides would be interchanged.) The increase in period (decrease in frequency) for this deep-fluid example is indicated by shifts in both the maximum am-plitude times and in the times for return to zero amplitude. Note especially that neither the two maxima nor the two zeros exactly coincide between the spike and bubble sides. In particular, at the end

of a half period, the interface is not perfectly flit. This result, together with the experimental results

of Fultz,'2 indicates a lack of perfect periodicity

in the fluid motion. Unfortunately, the computer

runs are too time consuming to be extended through

many complete periods, by which the statistics of

the oscillations could be examined.

Our principal interest was in verilyixg the shift of resonance frequency with amplitude and,

par-ticularly, to determine the depth-to-wavelength ratio at which the frequency shift reverses. it is apparent from Fig. 11, however, that from a hail-period

cal-culation four slightly different frequencies can be

calculated. Since the spike is relatively weliresolved,

especially in the shallow-fluid calculations, we in-vestigatd only the frequencies assóciated with the

spike mtion. To express our results, for comparison, in the n+tation of Fultz,'2 we define a dimensionless linear-tFeory amplitude

rn (TÄ/L)(/gD)1, (9)

03

FIG. 12. The

val-ues of , and w ea functions of h. Dot datum pointa are

based upon first manimum of spike motion; cross

da-turn points are

based upon first zero of spike

mo-tion. Plus datum

point shows

corn-bined spike and bubble result,

cor-rected to Fuitz'

definition of e.

SARY

-Comparison of the iesults of these calculations with those of experiments and analysis serves to demonstrate once again the validity of the marker and cell computing techñique for free-surface in compressible fluid flows. It also shows how cal-culations which include the full nonlinear features

of the Navier-Stokes equations in several space dimensions serve as controlled experiments, whose

detailed output can guide the develópment of

analysis as well as the planning of experimental

programs. -.

and a dimensionless measure of depth

h = vD/L. (10)

In terms of these, the frequency of oscillation is expressed in the form

w = w, + e2(kg)+wz,

so that w2, which is a dimensionless function of h,

is a measure of the deviation from the basic resonance

frequency o, of Eq. (8). The experimentaP2 and

theoretical" results for w2 (h) are reproduced in Fig.

12. Also shown are the calculation results, based

on the spike motion for the first half cycle. For the

deep-fluid example, the deviations from both the

theory and the experiments are relatively large.

For that case, however, we examined the effects

through a full period, by combining bóth the spike

and bubble frequency shifts. The result is much

nearer to the experimental curve.

Even the remaining discrepancy can be removed. Fultz'2 used a slightly different definition of e from

that of our Eq. (9). (Our definition follows that of

Tadjbakhsh and Keller.") Fuitz did not have

avail-able the necessary data for obtaining exactly our

value of ; the amplitude he used for computing it was obtained from the observed amplitude of the

oscillations, rather than from the linear-theory

am-- plitude. If we use this in our calculations for the full period of the deep-fluid example, we obtain

= 0.29, which compares favorably with the value from Fultz,w2 = 0.31.

The calculations of spike motion were ex'unined

further to clarify this relationship between actual

amplitude P and linear-theory amplitude Pò the latter being the amplitude of Eq. (8). Using the definition of e in Eq. (9), we found that

P/P0=1+,,e,

where the dimensionless j varies only slightly with h in a manner exhibited in Fig. 12.

Thus, for example, we have determined the bub-ble radius of curvature for both low- and

high-viscosity RayleighTaylor instability. No analysis known to the authors has succeeded in predicting

these results, and no previous experiments have been available for accumte comparison.

The calculations go much further, however, in

furnishing data for comparison with analysis. They show the nonlinear transition process in detail,

carrying several examples from the initial linear phase well into the final spike-and-bubble phase. They show the strongly modifying effects of

vis-cosity, both in speed of motion and in shape of

interface. In every case, the quantitative accuracy

THE PHYSICS OF FLUIDS

of the features presented in the figures has been demonstrated by comparisons with theoryor

experi-ment wherever available.

In the case of the nonlinear oscillatory motion, the frequency-reversal phenomenon is confirmed, and considerable additional data are made avail-able to show the time history of the motion for a particular type of initial impulse to the surface.

ACKNOWLEDGMENTS

We wish to acknowledge several valuable

con-tributions to this study by John P. Shannon. This work was performed under the aùspices of

the United States Atomic Energy Comini4sion.

L INTRODUCTION

IN

analytic solutions of the NavierStrokestextbooks and encyclopedias,, time-dependent equa-tions are compiled only in a sporadic way for special problems. The most comprehensive collection _can be found in Lamb's classical volume on

hydro-dynamics,' which may be supplemented by more recent books like those of Schlichting,2 Rosenhead,3

Landau, and Lifshitz,4 Dorfman,5 and Langl.ois.° 1 H. Lamb, Hydrodynamics (Dover Publications, Inc., New

York, 1945) 6th ed., Chap. Xl.

'H. Schlichting, Boundary Layer Theory (McGraw-Bill

Book Company, Inc., New York, 1960), 4th ed., Chap. V.

'L. Rosenhead, Laminar Boundary Layers (Clarendon Pressi Oxford, England, 1963), Chap. III.

EL.. D. Landau and E. M. Lifshitz, Fluid Mechmics (Peramon Press, Inc., New York, 1959), Chap. II.

'L. A. Dorfman, Hydrodynamic Resistance and the Heat

Loaa of Rotating Solids (Oliver and Boyd, Edinburgh, 1963).

S% E. Langlois, Slow Viscous Flow (The Macmillan Company, New York, 1964), Chap. IV.

VOLUME 9. NUMBER 5

Birth and Decay of Vortices

Hass J. Lu'r aND ERNST W. SCHWWERSKI U. .5. Naval Weapons Laboratory, Dahlqren, Virginia

(Received 17 May 1965; final manuscript received 29 December 1965)

Complete sets of time-dependent solutions for Stokes' slow-motion equations are constructed by means of a generalized separation technique, which was recently introduced and applied_{by the} authors to steady-state flows. Under certain boundary and regularity conditions, such integrals of slow motion can be used to derive solutions of the quasi-linear Navier-Stokes_{equations in an exac,,}

linear manner. The classes of separable solutions, which are presented in this paper, reveal significant

mathematical properties of laminar motions concerning their uniqueness and regularity. From the physical point of view, they render valuable information on the mechanism of growing and dying yortices. For instance, it can be shown that a decaying discontinuity line is _{a process which involves} the gradual obliteration of a discrete spectrum of vortices in such a way that little vortices with large damping coefficients dominate first, and that bigger vortices of fewer number replace them with advancing time. Furthermore, the development and decay of singular vortices are investigated.

MIAY 1966

Most of the familiar integrals of time-dependent flows share the common property that they fulfill

simultaneously Stokes' linear equations. (An ex-ception is Taylor's special solution of the

Navier-Stokes equations,7 which describes a symmetric

net-work of decaying vortices of equal strengths.) _As

verified by inspection, those solutions do not show

any development of secondary vortices. Instances for rotating flows are HamelOseen's solution of a decaying potential vortex' and the accelerated

mo-tion of a circular cylinder.3 Stokes' first and_second

problems of suddenly accelerated and oscillating plates2 and a decaying straight discontinuity _line'

are examples of parallel flows. For time-dependent G. I. Taylor, Aeronautical Research Communication (1923); see, also, G. K. Batchelor, The Scientific Papers of Sir Geoffrey Ingram Taylor (Cambridge University _Press,