The full Navier-Stokes equations are solved numerically in cylindrical coordinates in order to investigate the splash of a liquid drop onto a flat plate, into a shallow pool, or into a deep pool. Solution is accomplished with the Marker-and-Cell technique using a high-speed computer. Results include data on pressures, velocities, oscillations, droplet rupture, and the effects of compressibility. They also show how the technique can be applied to a wide variety of other complicated fluid flow problems involving the transient behavior of a free surface.
Received 17 May 1967
The remarkable appearance of a splashing liquid drop has been observed for centuries, but detailed studies of the process were not possible until the development of high-speed photographic techniques. The first extensive description was given by Worthington [1]. More recently the phenomenon has been investigated for a variety of specialized purposes; Edgerton and Killian [2] used it to illustrate their photography techniques; Gillespie and Rideal [8] were concerned with the bouncing of drops from aerosol filters; Sahay [4] examined the breaking of water drops falling into oil; and Hobbs and Kezweeny [5] studied the breakup of the rebounding splash jet.
These previous investigations are principally experimental. Full theoretical studies have been precluded by the great difficulty in solving the appropriate equations. Our approach has therefore been numerical, using the Marker-and-Cell technique [6] [7] [8] to obtain high-speed computer solutions of the full Navier-Stokes equations for a viscous, incompressible fluid. The purpose is to determine the flow dynamics during the splash, and thereby to explain some of the experimental results that have been published with little or no correlative interpretation.
The splashing properties of a liquid drop of density ρ
are determined by several dimensionless parameters
that are functions of the drop radius R,
the collision speed u0,
the depth of the splash pool D,
the acceleration of gravity g,
the kinematic viscosity coefficient ν,
and the surface tension coefficient T.
(We assume that the drop is spherical and that the fluid in the drop is
the same as that in the pool.) From these we form the four
dimensionless parameters,
R/D,
(gR)½/ua, u0,
u0R/ν and
T/ρRu02,
which characterize the flow.
If the drop had fallen with negligible air resistance from a height h, then the second dimensionless parameter would equal (R/2h)½. Actually, air resistance retards the drop, ultimately to a terminal speed that depends upon drop size, but does not exceed about 10 m/sec. The air tends also to deform the drop, and it may even oscillate from prolate to oblate during fall, but this effect is negligible for drop radii less than about 0.15 cm [9]. Surface tension, which serves to maintain sphericity, also enters strongly into some aspects of the collision dynamics, but our calculations neglect these effects. If T/ρRu02 is sufficiently small (a condition easily achieved with water drops) then the principal effects of surface tension occur at the edges of the thin lateral sheet jets, which are of little concern to the principal purposes of our calculations. The inclusion of surface tension in Marker-and-Cell calculations, however, has been described by Daly and Pracht [10] for some similar investigations.
The Marker-and-Cell technique for high-speed computer has been described in several publications [6] [7] [8], so that only a brief description is given here.
The full (nonlinear) Navier-Stokes equations in cylindrical coordinates (with axial symmetry) are approximated by finite-difference expressions corresponding to a mesh of computational cells over the domain. Using the finite difference incompressibility condition, we obtain a Poisson equation for the pressure. The initial-value problems are solved by advancing the configuration through a set of finite time steps, or cycles. Each cycle consists of four steps:
The results, like those of an experiment, are recorded periodically. These include the configuration of marker particles, the velocities, pressures, and any desired functionals of those quantities. The illustrations of particle configurations and pressure surfaces in this paper were processed directly from the computer by the Stromberg-Carlson SC-4020 microfilm recorder and have not been retouched or otherwise altered. Other functions have been graphed from printed data. The calculations themselves were performed by a code in Fortran IV on the CDC-6600 computer.
The numerical stability of these calculations has been discussed by Daly and Pracht [10], extending the concepts described in the previous Marker-and-Cell method references, and illustrating an example of the powerful non-linear stability theory of Hirt [11].
The accuracy has been demonstrated by the usual tests of varying cell size and computational cycle interval, and insisting that these be small enough for negligible difference in the results.
No novel methodological features have been introduced beyond those previously described. The use of cylindrical coordinates is a straightforward extension of the technique as developed for Cartesian coordinates. (We now believe that even in Cartesian coordinates it is useful to base the difference equations on the viscous term −ν∇×∇×u rather than ν∇2u, as previously described. The result is considerable simplification, rather than increased accuracy.) The Marker-and-Cell technique has now been generalized to incorporate the correct free surface stress conditions, but these are not well demonstrated in the present studies, in which the effects of viscosity are negligible. Recent calculations by J P Shannon on the Coanda effect and by Hirt and Shannon on hydraulic jumps demonstrate the tremendous importance of the viscous contributions to the free surface stress for flow problems with appreciable viscosity.
It should be emphasized that the Marker-and-Cell technique is equally applicable to confined-flow problems, and to circumstances in which there are two or more fluids of different densities [8]. Even in those examples, the technique appears to have advantages over others that use the vorticity and stream function as primary variables [12].
In all calculations described in this paper, the effects of surface tension and viscosity have been neglected.
Accordingly, the results are completely determined by the parameters R/D and (gR)½/u0. The splash onto a flat plate is characterized by R/D=∞.
Figure 1 shows a typical sequence of results. The magnitude of (gR)½/u0 is zero. As observed experimentally, such a drop does not splash upwards, but instead forms a lateral sheet jet.
The calculation shows that the jet moves with constant speed, 1.6 times greater than the initial impact speed of the droplet. Actually the results of Worthington [1] show that the lateral jet speed decreases with time, presumably because of the surface tension. Measurements from Worthington's figures (Ref 1, p. 137) show, nevertheless, rather close agreement with the computer results until the time t=12, as shown in Fig 2. Thereafter the experimental jet tip slows down and eventually is pulled back inwards by the surface tension.
The calculation shows a second outward moving wave with speed u0; this is the shoulder that shows well in the latest-time frame of Fig 1. The flow behind this shoulder is like that which occurs in an impinging jet, for which the Bernoulli theorem indicates a lateral velocity equal to u0 in steady state. The novel feature here, therefore, is the thin sheet of much faster fluid that is initiated during the first stages of impact, when steady-state theories are not applicable. With its relatively great energy density, the fluid in this sheet jet is capable of considerably enhanced erosive or cutting power, an important consideration in agricultural engineering [13].
The calculation also represents the collision of two drops that impact along the line of centers. The lateral sheet jet speed would then be 0.8 times the initial relative collision speed.
When the pool depth is comparable to the drop radius, then the splash properties vary rapidly with R/D. To emphasize these we first examined a number of cases for which (gR)½/u0 is negligible small. For all calculations in this section, R=10.0 and u0= 1.0, thereby establishing the units of time, pressure-to-density ratio, etc.
The striking feature is that for R/D large but finite, the splash is distinctly different from that for R/D=∞ . Even a very thin film of fluid in the pool is enough to interact appreciably with the lateral sheet jet to produce an upward motion. This can be seen in the photograph by Edgerton reproduced in Fig 3, and in the shallow-pool calculation shown in Fig 4, for which R/D=5.0. The contrast with Fig 1 is considerable.
One feature of disagreement between Fig 3 and Fig 4 is seen in the experimental tendency to form droplets at the edge of the splash sheet. This is precluded in the calculations by the ommision of surface tension and by the restriction of axial symmetry. (Three-dimensional calculations cannot be resolved sufficiently by present-day computers to make the full calculation practical. The Marker-and-Cell method, however, is completely applicable to three-dimensional investigations without requiring any basic modifications.)
Figure 5 shows another calculation; in this case R/D=1.0. Even for this great a pool depth, the tendency to empty the region of impact is quite marked.
Many of the properties of the shallow-pool splash are shown by the sequences of fluid configurations. Figure 6, for example, illustrates several of the features that can be measured. Notice that in the absence of gravity, the lateral splash sheet grows in height and radius with uniform speed after an initial establishment period. This is true for all of the calculations with R/D in the interval from 1.0 to 5.0. The results can be expressed for the rate of radius increase, ur, and of edge height rise, ue, as follows:
ur/u0 = 0.77 − 0.20 (D/R)2
ue/u0 = 0.36 − 0.05 (D/R)2
The results for a shallow pool are in agreement with those that can be measured from the photographs of Edgerton and Killian [2], from which we obtain ur/u0=0.8±0.1 and ue/u0=0.4±0.1.
In addition, there are properties of the flow, such as pressure, that cannot be measured directly from the fluid configurations. Figure 7 shows pressure surface plots for the same problem as illustrated in Fig 5, and for the same times as the last two frames shown in that figure. In Fig 7, the vertical axis in the distance is the ordinate of pressure, whose base rests at the intersection of the cylindrical axis with the bottom of the pool. Radiating from that base are the cylindrical axis to the right and the bottom of the pool to the left. Each plot is scaled relative to the maximum fluid pressure at that instant; this variation with time is shown in Fig 8.
Another plot of pressures is given in Fig 9, which shows pressure-to-density ratio along the bottom of the pool at various times for a calculation in which R/D=2.0. From data of this sort impulse distribution on the bottom can be calculated in whatever form one desires.
For this same calculation, the radial velocity along the bottom is plotted for several times in Fig 10. These results can be used to estimate shearing stress and drag erosion, provided that the boundary layer is sufficiently thin.
For a deep pool, the crucial splash parameter is (gR)½/u0. With present computers, calculations can be performed only for relatively large values of the parameter; for low values (high impact speed) the early-stage crater size is too large to be resolved.
The three examples chosen for illustration have R=5.0 and g=0.1, and initial impact speeds such that (gR)½/u0=0.707 (Fig 11), 0.354 (Fig 12), and 0.177 (Fig 13).
These correspond, respectively, to fall distances of R, 4R, and 16R. At the extremes, they demonstrate quite different types of droplet distortion. In Fig 11, the original droplet fluid maintains its integrity, whereas in Fig 13 the droplet is broken into two parts, one remaining below the surface in a fairly compact volume, the other spread out as a film on the surface.
Sahay [4] has observed this distinction experimentally, and has determined the critical falling distance as a function of drop size for water falling into viscous oils of lower density than water. For shorter falls than the critical, the water drops retain their integrity during splash; otherwise they break into two. Sahay shows that the critical fall distance is inversely proportional to the droplet mass. In the absence of viscosity or surface tension, the scaling parameter shows that the critical height is proportional to the cube root of the droplet mass, while if viscous forces are dominant it is inversely proportional to the two-thirds power of droplet mass. Thus, the results of Sahay indicate a strong effect of viscosity, in contrast to the results presented here. Nevertheless, we believe that the basic mechanism responsible for the rupturing of splashing drops is as we present it here, and that viscosity serves only to modify the quantitative aspects by giving the droplet an apparently smaller impact speed.
In all three illustrated examples, the initial stages are quite similar. The droplet material is deposited onto the bottom and sides of a crater, the bulk of the material lying near the bottom. As the cavity sides collapse, the bottom simultaneously tends to rise back towards the original surface level. Thus, two processes compete. For a low-speed impact (Fig 11) the crater is shallow and the droplet rebounds to the surface before the sides collapse. For a high-speed impact (Fig 13), the crater is much deeper and the sides collapse earlier, trapping much of the original drop material well below the surface. In this latter case, the material deposited on the sides of the crater is carried up in the rebound splash, isolating two different parts of the original drop fluid from each other.
If the pool fluid were the less dense and the two were immiscible, then each part of the droplet would coalesce and the two would drift slowly downwards. In our calculation, the larger of the two parts would be the lower, so that they would not tend to reunite. In the absence of any comments in this regard by Sahay [4], we expect that this explanation pertains directly to his observations.
Worthington [1] implies that the drop material is all contained in the rebounding central column. The experimental evidence is based upon photographs of milk drops falling into water and of droplets coated with a thin layer of soot. Our calculations show that part of the original droplet fluid does indeed rise as a coating to the rebound jet, thereby agreeing with Worthington's observations. But much of the droplet fluid can be trapped in a pocket below the surface, and we suppose that this was not visible in Worthington's experiments because of the refractive properties of the water.
Figure 14 illustrates the motion of the top and bottom of the droplet for the calculations shown in Fig 11 and Fig 12. The oscillatory motions are related to those of the Cauchy-Poisson problem discussed by Lamb [14] and Stoker [15], 5 but an analytical treatment of low-impact-speed splash problems appears to be much more difficult. Estimation of the Cauchy-Poisson period of oscillation for the examples of Fig 14 gives results that are approximately one-third of those indicated at the latest times in the figure.
The heaving oscillations of a freely floating body are also related. John [16] has discussed this for some simplified examples, obtaining results that qualitatively resemble those in Fig 14. Again, a really pertinent analysis for the splash problem appears intractable.
The impact of certain plastic or metallic spheres may produce melting, in which case the splash resembles the incompressible processes already discussed. As the impact-speed increases (approaching the order of 5 mm/sec) the effects of compressibility become progressively more important; the material vaporizes during impact and subsequently behaves much like a gas. An example to show this extreme is given in Fig 15. The material is a polytropic gas with γ=5/3. The radius of 10.0 and impact speed of 1.0 are the same as in Fig 4 and Fig 5, but the contrast with those figures shows several important effects of compressibility.
At first there is almost no splash. Shocks are generated in both the drop and the pool, but their lateral motion in the pool is relatively slow so that very little expansion occurs in the surface. As soon as the shock reaches the top of the droplet, there is a strong reexpansion that violently expulses the material back out in the direction from which it entered. The outward moving speed of the droplet top is very close to the initial impact speed.
This set of calculations was performed by A A Amsden, using the Particle-in-Cell technique for high-speed flows [17] [18].
We are grateful to J E Welsh for the Marker-in-Cell computer code, to A A Amsden for the compressible flow calculations, and to H E Edgerton for the splashing drop photograph. This work was performed under the auspices of the United States Atomic Energy Commission.
[1] A M Worthington, A Study of Splashes, MacMillan, (1963).
[2] H E Edgerton, J R Killian Jr, Flash Charles T Branford, Boston, (1954).
[3] T Gillespie, E Rideal, J Colloid Sci. 10, 281 (1955).
[4] B K Sahay, Indian J Phys, 18, 306 (1944).
[5] P V Hobbs, A J Kezweeny, Science 155, 1112 (1967).
[6] F H Harlow, J E Welch, Numerical Calculation of Time-Dependent Viscous Incompressible Flow of Fluid with Free Surface, Phys Fluids 8,2182 (1965).
[7] J E Welch, F H Harlow, J P Shannon, B J Daly, The MAC Method, Los Alamos Report, LA-3425, (1967)
[8] B J Daly, Phys Fluids 10,297 (1967).
[9] D M A Jones, J Meteorol 16 504 (1959).
[10] B J Daly, W E Pracht, Phys Fluids, submitted for publication.
[11] C W Hirt, J Computational Phys, submitted for publication.
[12] J E Fromm, F H Harlow, Phys Fluids 6, 975 (1963).
[13] A A Amsden, The Particle-in-Cell Method for the Calculation of the Dynamics of Compressible Fluids, LASL Report NLA-3466 (1966) .
[14] H Lamb Hydrodynamics, Dover, pp 384,430 (1945).
[15] J J Stoker, Water Waves, Interscience, p 156 (1957).
[16] F John, Comm on Pure and Appl Math, 2, 13 (1949).
[17] F H Harlow, Proceedings of the Symposia in Applied Mathematics, American Math Soc, Providence, Vol 15, p269 (1963).
[18] A A Amsden, The Particle-in-Cell Method for the Calculation of the Dynamics of Compressible Fluids, LASL Report LA-3466 (1966).