Received 8 August 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 out 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, non-linear Navier Stokes equations for an incompressible fluid.
WHEN a box of fluid has been subjected to an impulsive force on its free surface, the nature 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 case, 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 [1] 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 dimensions, 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 computations are based, by direct comparisons of the results with physical experiments. Second, they give information about the processes for 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 complement the work of the laboratory experimenter and of the analytical theoretician, but certainly do not replace either.
All calculations reported in this paper were performed on the IBM 7030 (Stretch) Computer. The figures showing particle configurations, velocity vectors, or pressure contours were processed directly from the computer by the Stromberg-Carlson SC-4020 Microfilm Recorder and are not retouched or otherwise altered.
The marker and 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 δx by δy. This enables the fluid to be represented by a finite set of values of the appropriate. field variables. In the marker and cell method, these variables are primarily the pressure and the components of velocity. 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 showing 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 addition, the time derivatives are similarly approximated, so that the calculation proceeds from the initial conditions through a sequence of cycles, each of finite duration δt. The finite-difference approximations are derived from the Navier-Stokes equations written in the following conservative form:
The velocity components u and v are parallel, respectively, to the x (horizontal) and y (vertical) axes. The ratio of pressure to (constant) density is designated by φ, ν is the kinematic viscosity coefficient; and gx, gy, are the components of body acceleration (gravity). The only approximations involved in the solution 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:
These four steps complete the basic parts of a cycle of time advancement. Actually, there are many minor aspects, presented previously [1], 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 computational cells be sufficiently fine to resolve the essential details of the flow field. Also, the time step per cycle (δt) must be small enough to ensure that field-variable changes per cycle are almost everywhere much less than the values of the variables. (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 δt is sufficiently small if the projected viewing appears 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 consideratons of calculational 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 possibility of numerical instabilities, of which the investigator must always be wary. For the marker and cell calculations, this consideration 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 considerable accuracy when ν=0. The 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 ν=0, the irregular results may appear turbulent, 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 ν 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 are 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. The results shown in this paper give only a small sample.
With gravity pointing upward, out of the fluid, the free-surface irregularities are unstable. Thus, after an initial sinusoidal pressure impulse on the surface, the fluid continues to fall upward with increasing 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 Birkhoff [5] [6], 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 analytical techniques are readily available, Such modifying features as surface tension, viscosity, and continuous layering have also been considered.
Experiments show, however, that when the perturbation 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 presented by Birkhoff [6]), by series expansions of limited validity range (see Emmons, Chang, and Watson [8]), or by numerical techniques. Some of the numerical methods for incompressible fluids have been summarized by Birkhoff [5] [6], but none have been extensive enough for useful analysis.
The present technique,. while applicable to only one fluid with a free surface, enables computer calculations 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 scaled in such a way that the magnitude of the body acceleration is unity and the half wavelength of the 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/ρ≡φ, 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
φs=Aδt cos kx
in which k=π/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 the 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 1-5. For this example, A=5.0 (large impulse), D=8.0 (deep water), and ν= 0.01 (negligible viscous effect). Figure 1 shows the marker particle positions at a sequence of times, illustrating 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 local fluid speed, and whose direction is the local flow direction.
The lines in Fig 3 are isobars, which show the dramatic changes occurring in the pressures at early times. A peculiarity of the plot program slightly shifted the isobars 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(gξ/D)½] cos kx (1) in which ξ ≡ kD tanh (kD) (2) For early times, Z = −(Aξt/D) cos kx, (3) referred to as the initial impulse motion.
For the calculation described in Figs 1-3, Eq (1) applies 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 displacement histories of the spike and the bubble. The line splitting the two at early times shows the initial 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 curvature of the bubble. Birkhoff [5] [6] has discussed studies of several investigators concerning bubble motion, and jointly with Carter [6] [9] has developed an analytical prediction for steady-state bubble speed. The predicted result [6] is indicated by one of the horizontal lines in Fig 5. A similar result by Garabedian [10] 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 Birkhoff [5] suggests that his accelerations were considerably greater than reported, thereby making bubble-speed analysis uncertain.
The radius of curvature of the bubble in our calculation, measured in units of perturbation wavelength, 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 experimental 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 described 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 somewhat 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 freefall stage, as expected.
For these calculations with negligible viscosity, the accuracy of the results has been assessed principally by comparisons with theory and experiments at the earliest or latest phases of the motion. For the intermediate transitional phase, there is.very little comparative material. One exception is furnished by the calculations of Emmons, Chang, and Watson [8], whose theoretical treatment included effects to third order in the amplitude. Graphs of their interface configuration show clearly the beginning 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 their neglect of higher-order effects contributes to an inaccuracy which, in comparison with the computer results, shows a slightly narrower spike and broader bubble.
It is well known that the effects of viscosity and surface tension on Rayleigh-Taylor 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 competition between two energy sources. The gravitational body force accelerates the flow as before. The viscous dissipative forces, however, tend to slow down the flow in a manner which can be described, for the extreme case in which g=0, by the surface behavior
Z = − (A cos kx/kνψ) [1 − exp(ψk2νt)] (4)
The dimensionless number ψ is the nonzero root of the equation
ψ2 + 4 ψ = 4[(1+ψ)½ − 1] (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 viscous damping must dominate over the gravitational acceleration.
More generally, the full eigenvalue equation for the exponential growth rate coefficient can be written
1 − ξ2 − 4 αξ + 4α2 {[1+(#958;/α)] ½ − 1} = 0 (6)
in, which ξ and α are dimensionless measures of the time coefficient and of the viscosity, respectively,
ξ2 ≡ n2 g-1k-1, α2 ≡ k3 ν 2 g-1 (7)
[The amplitude is proportional to a weighted sum of terms in exp (nt).] 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 α2 ≪ 1 ie, if ν2≪ gk-3. Under that circumstance, the appropriate solutions of Eq (6) are ξ= 1−2α+2α3/2 and ξ=−1−2α−2iα3/2, with fairly good accuracy for α≤0.1. For large values of α, the real root approaches (2α)-1, while the complex one is proportional to α, 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 illustrated in Figs 1-5, except that the coefficient of viscosity is ν=2.0. Correspondingly, α=1.058, and the real root of Eq (6) is ξ=0.377. [In contrast, the low-viscosity calculations, with ν=0.01, have α= 0.0053 with only a 1.1% effect on the roots of Eq (6).]
A comparison between Fig 1 and Fig 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 differences shown by the comparison. Viscosity increases the width of the spike and narrows the bubble, so that the basic harmonic of the perturbation persists longer.
Comparison with Fig 2 shows that there also are qualitative differences in the appearance of the field of velocities. Strong motions of convergence into the spike characterize the low-viscosity case, while the large-viscosity effect is to make the velocities near the spike center more nearly vertical.
A quantitative comparison of the viscous effect is visible in
Fig 4 and Fig 8.
(Note the necessary difference in time scale between the two.)
In Fig 4 the theoretical initial impulse solution
lies between the calculational results for bubble and spike, whereas
Fig 8 shows that the strong effect of viscosity
is sufficient to retard both spike and bubble motion to well below the
initial-impulse curve. Also shown in Fig 8 is a
plot of Eq (4); the close agreement with the mean 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
uB(2Lg)-½=0.184,
while its radius of curvature is given by R/2L=0.313. These, which are
nearly the final asymptotic values, may be compared 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 of the process. Thus, this same large-viscosity run was repeated
with a much smaller initial impulse, A=0.5. Fig 9 shows the
full history of the configuration development, 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 to the acceleration g=1.0, while the
bubble is characterized by
UB(2Lg)-½=0.181
and
R/2L=0.303.
(These values confirm the statement that the calculation shown in
Fig 8 has nearly reached its steady state.)
Theoretical calculations [11] and confirming experiments [12] 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 to wavelength ratio is sufficiently small, there is an increase in frequency.
To investigate this effect, we changed the previously-described calculations only insofar as directing the body acceleration into the fluid. All scaling and dimensions remained unchanged; i.e., L=4.8 (half wavelength), g=1.0, ν=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.
The results of linear analysis predict a surface motion
Z = −A(ξ/gD)½ sin [t(gξ/D)½] cos kx (8)
which also is plotted for comparison in Fig 11. 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 so 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 amplitude 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 flat. This result, together with the experimental results of Fultz [12] 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 verifying the shift of resonance frequency with amplitude and, particularly, to determine the depth-to-wavelength ratio at which the frequency shift reverses. It is apparent from Fig 11, however, that from a half-period calculation; four slightly different frequencies can be calculated. Since the spike is relatively well resolved, especially in the shallow-fluid calculations, we investigated only the frequencies associated with the spike motion. To express our results, for comparison, in the notation of Fultz [12], we define a dimensionless linear-theory amplitude
ε ≡ (πA/L)( ξ/gD) ½ (9)
and a dimensionless measure of depth
h= πD/L (10)
In terms of these, the frequency of oscillation is expressed in the form
ω = ω0 + ½ε2(kg)½ ω2
so that ω2, which is a dimensionless function of h, is a measure of the deviation from the basic resonance frequency ω0 of Eq (8). The experimental [12] theoretical [11] results for ω2(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 both the spike and bubble frequency shifts. The result is much nearer to the experimental curve.
Even the remaining discrepancy can be removed. Fultz [12] used a slightly different definition of ε from that of our Eq (9). (Our definition follows that of Tadjbakhsh and Keller [11]). Fultz did not have available 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 amplitude. If we use this in our calculations for the full period of the deep-fluid example, we obtain ω2=−0.29, which compares favorably with the value from Fultz, ω2=−0.31.
The calculations of spike motion were examined further to clarify this relationship between actual amplitude P and linear-theory amplitude P0, the latter being the amplitude of Eq (8). Using the definition of ε in Eq (9), we found that
P/P0 = 1 + ηε
where the dimensionless η varies only slightly with h in a manner exhibited in Fig [12].
Comparison of the results of these calculations with those of experiments and analysis serves to demonstrate once again the validity of the marker and cell computing technique for free-surface incompressible fluid flows. It also shows how calculations 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 development of analysis as well as the planning of experimental programs.
Thus, for example, we have determined the bubble radius of curvature for both low- and high-viscosity Rayleigh-Taylor instability. No analysis known to the authors has succeeded in predicting these results, and no previous experiments have been available for accurate 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 viscosity, both in speed of motion and in shape of interface. In every case, the quantitative accuracy of the features presented in the figures has been demonstrated by comparisons with theory or experiment wherever available.
In the case of the nonlinear oscillatory motion, the frequency-reversal phenomenon is confirmed, and considerable additional data are made available to show the time history of the motion for a particular type of initial impulse to the surface.
We wish to acknowledge several valuable contributions to this study by John P. Shannon.
This work was peformed under the auspices of the United States Atomic Energy Commission.
[1] F H Harlow, J E Welch, Phys Fluids 8,2182(1965).
[2] F H Harlow, J P Shannon, J E Welch, Science 149, 1092 (1965).
[3] G I Taylor, Proc Roy Soc, A201, 192 (1950).
[4] D J Lewis, Proc Roy Soc, A202, 81 (1950).
[5] G Birkhoff, LASL Report LA-1862 (1955).
[6] G Birkhoff, LASL Report LA-1927 (1956).
[7] S Chandrasekhar, Hydrodynamics and Hydromagnetic Stability, OUP, 1961.
[8] H W Emmons, C T Chang, B C Watson, J Fluid Mech 7,177 (1960).
[9] G Birkhoff, D Carter, J Math and Mech, 6,769 (1957).
[10] P R Garabedian, Proc Roy Soc, A241,423 (1957).
[11] I Tadjbakhsh, J B Keller, J Fluid Mech, 8,442 (1960).
[12] D Fultz, J Fluid Mech, 13,193 (1962).