Received 16 April 1965, Final manuscript 3 September 1965
A new technique is described for the numerical investigation of the time-dependent flow of an incompressible fluid, the boundary of which is partially confined and partially free. The full Navier-Stokes equations are written in finite-difference form, and the solution is accomplished by finite-time-step advancement. The primary dependent variables are the pressure and the velocity components. Also used is a set of marker particles which move with the fluid. The technique is called the marker and cell method. Some examples of the application of this method are presented. All non-linear effects are completely included, and the transient aspects can be computed for as much elapsed time as desired.
THIS paper describes in detail a new technique for the numerical solution of problems in the dynamics of an incompressible fluid with a free surface. The method has been developed for use on a high-speed electronic computer and would be impractical for hand-solution purposes. It is here designated the marker and cell method.
Specifically, the problems to which the method applies involve the time dependent motion of a viscous, incompressible fluid in two-dimensional Cartesian coordiuates. The fluid may be bounded in part by the walls of an irregular box or by lines of reflective symmetry. A uniform body force (gravity) may act on the fluid, and a prescribed pressure may be applied to the surface, variable in both position and time. The investigator supplies to the computer the initial and boundary conditions for a specific study of interest, and the computer then develops the solution ihrough a succesion of closely spaced time increments for as long as the results continue to be of interest. The analogy to a physical experiment is therefore quite close.
The solution technique makes use of finite-difference approximations applied to the full Navier-Stokes equations. The primary dependent variables are the pressure and the two components of velocity. Neither the vorticity nor the stream function enters explicitly into the analysis.
The finite-differences apply to both space and time variations. For the former the region is divided into numerous small rectangular zones or cells and the field variables in each are characterized by single average values. For the time variations the changes are represented by a succession of field variable values separated by small increments of time. If both the finite space and time differences are small enough, the results will be sufficiently close to continuous while large differences may destroy essential resolution to the point of producing nonsense.
In recent years there have been many techniques developed for the numerical solution of complicated problems in fluid dynamics. Spectacular results came first for compressible (high-speed) flows, and a large literature has been produced on the subject [1]. Incompressible flows are now also being studied numerically with considerable success, and new methods rapidly are becoming available.
One of the principal applications of incompressible flow calculations has been to the problem of weather analysis and prediction [2]. Problems of natural convection have also been extensively investigated [3] [4] [5] [6]. A useful survey of previous work has been given by Pearson [7], who has also proposed some computing techniques. A group at Wisconsin [8] is now undertaking the study of a variety of additional numerical approaches to the solution of the Navier-Stokes equations. The published reports of these more recent investigations contain bibliographies from which can be traced much of the earlier work.
At Los Alamos the efforts in computing incompressible fluid dynamics have been directed primarily to nonsteady viscous flows [9]. The principal published applications have been to the von Karman vortex street behind a cylinder [10] [11]. As in almost every other incompressible-flow technique, the primary dependent variables used in the calculations are the vorticity and the stream function.
For the case of free surface flows, we have found in the marker and cell technique that the primary physical variables, velocity and pressure, have several advantages over the stream function and vorticity. The free surface boundary condition of vanishing stress, or of prescribed normal stress, becomes more natural to apply. In addition, the physical significance of any new modification in technique is more readily visualized.
The detailed derivation of the finite difference equations is based upon the following sequence of events by which the configuration is advanced from one time step to the next.
This, then, completes the advancement of configuration to the end of the new cycle. Actually, several crucial points have been omitted in this brief summary; these are more easily discussed in the detailed description to follow.
One point should be emphasized. The marker particles introduced into this incompressible flow calculation are only for the purpose of indicating fluid configuration. They show which cells contain fluid and especially which cells lie along the free surface. They also serve as a flow visualization coordinate system whereby fluid element trajectories and relative positions can be observed. They do not participate in the calculation as do the Particle-in-Cell-method particles in the Los Alamos computer program for compressible flows (see Ref 1, pp. 319- 343).
The starting point for the derivation is the Xavier-Stokes equations for viscous incompressible flow in two-dimensional Cartesian coordinates. We use a differential form which is completely conservative of mass and momentum,
The x and y coordinate axes are horizontal and vertical, respectively, with the origin located at the lower left of the computing region. The corresponding components of the velocity are u and v, while φ is the ratio of the pressure to the (constant) density. (For brevity, we usually refer to φ simply as pressure.) The kinematic viscosity coefficient is designated by ν. The body force is designated by the constant accelaration components gx, and gy.
The finite-difference approximation to the equations corresponds to an Eulerian mesh of cells covering the computation region, as sketched in Fig. 1, For each cell the local field variable averages are centered as shown in Fig. 2. The cells are numbered by indices i and j which count cell center positions along the horizontal and vertical directions, respectively. Cell boundary positions are labeled with half-integer values of the indices. The dimensions of the rectangular cells are δx and δy.
In addition to the space index subscripts, we use a superscript to number the time cycle. For example,
designates the horizontal velocity at the time t = (n + l)δt, in which δt is the time increment per cycle. When the cycle number superscript is omitted, it is assumed that its value is n, corresponding to the time t = nδt. (The time cycle advancement to be discussed is from n to n + 1, so that the omission of a superscript refers to the value of the quantity at the beginning of the cycle.)
With these definiiions the finite-difference approximations to Eq. (2) and (3) can be written in the following form:
It may be noted that some of the velocity values are not centered at points indicated in the mesh diagram. Whenever this occurs an average of adjacent values is implied. Representative examples are
Define
Then the finite-difference analog of Eq (1) is Di,j=0.
From Eqs (4) and (5) we may obtain
an equation which is of fundamental importance to the derivation. We have used Qi,j, as an abbreviation for the following:
The procedure for determining pressures is based upon the requirement that Di,jn-1 vanish for every cell at the end of the time cycle. This leads to the equation for φi,j
where
Insertion of this into Eq (7) shows that Di,jn+1 ≡ 0 can result from calculation of the new velocities, provided that Eq (9) is accurately solved for the pressures. In principle, Qi,j could be used in Eq (9) instead of Ri,j, since the two differ by terms proportional to D. In practice, however, the use of Ri,j is desirable in that the solution of Eq (9) need not be nearly so accurately derived to keep the accumulation of compressional discrepancy to a sufficiently low level. Since Eq (9) is solved by an iterative procedure (see, for example, Ref 9), there is considerable economy in computer usage resulting from any process which decreases the required accuracy of solution. Tests have shown, however, that with a very stringent convergence requirement, the cumulative results of a calculation are independent of whether Qi,j or Ri,j is used in Eq (9).
In order that Di,jn+1 identically vanish for every cell, it is necessary also that Di,j vanish at all boundaries of the fluid. This requirement forms a useful basis for determining the finite difference analogy to the necessary boundary conditions.
The Eqs. (4), (5), (9), and (10) are the principal ones used in performing the calculations. First, the value of Ri,j is computed for every cell, using in Eqs (8) and (10) velocities available at the beginning of the cycle. Second, Eq (9) is solved for the value of φi,j for every cell. Third, the pressures so obtained are put into Eqs (4) and (5) ad the new velocity components arc computed.
If there were no free surface to the fluid, the essential parts of the calculation would then be complete, and after some bookkeeping processes the next cycle could commence. To include a free surface, one more step is essential. The new position of the surface must be calculated. This implies a scheme for keeping track of surface position. The one we have developed actually goes beyond this basic requirement; it supplies a coordinate system of marker particles whose trajectories follow the motions of elements throughout the fluid. Marker particles are initially placed in the cells containing fluid, and they subsequently are moved with local velocity. A linear interpolation is performed to calculate the velocity with which a particle should move. The interpolation weighting depends upon the distance of the particle from the nearest velocity points in the Eulerian mesh of cells.
A cell with no marker particles is considered to contain no fluid. A cell with marker particles, lying adjacent to an empty cell, is called a surface cell. All other cells with particles arc considered to be filled with fluid.
A rigid wall may be either of two types, no-slip or free-slip. The latter may be considered to represent a plane of symmetry, rather than a true wall, or if the fluid is viscous. it may represent a non-adhering (greased) surface. Walls are restricted in orientation so that they lie along the boundaries of the Eulerian calculational cells. Relaxation of this requirement could be accomplished only at the expense of considerable inerease in complication.
A vertical wall therefore passes through the horizontal-velocity mesh points, and the velocities at those points vanish at all times for either type of wall. A vertical wall does not pass through vertical-velocity mesh points, but the calculation makes use of the values of v at mesh points lying just outside of the wall. In the example shown in Fig 3, the calculations require a value of v', a typical vertical velocity beyond the wall. (In the following discussion, a prime refers to the exterior cell quantity.) For a no-slip wall the boundary condition is v' = -v, while for a free-slip wall it becomes v' = + v. Analogous boundary conditions are applied at a horizontal wall. In general, the tangential velocity reverses or remains unchanged according as the wall is no-slip or free-slip. The normal velocity reverses for a free-slip wall, wherever it is needed in Eq (4) or (5) or in the calculation of Ri,j. For a no-slip wall the normal velocities beyond the wall must be calculated in such a way as to ensure that D vanishes for the exterior cell.
Wall boundary conditions are also needed for the solution of the pressure equation. While it is not necessary to calculate the changes in normal velocity for points lying on the wall, the boundary conditions for φ must be consistent with the identical vanishing of that velocity. The derivations are accomplished using Eqs (4) and (5).
In the case of a free-slip wall, the boundary condition is easily derived. For a vertical wall refer to Eq (4). With reversal of all normal velocities and no change in all tangential velocities, every term vanishes except those with φ or g, leading to the boundary condition φ' = φ ± gxδx. The sign is + if the fluid lies to the left of the wall, - if it lies to the right. For a horizontal wall reference to Eq (5) shows a similar vanishing of terms, leading to φ' = φ ± gyδy where the sign is + if the fluid lies below the wall and - if the fluid lies above.
For a no-slip wall the velocitics beyond the wall must always be chosen in such a way that D' vanishes in order that D is prevented from diffusing into the fluid. This results in a somewhat more complicated set of pressure boundary conditions than for the free-slip wall. Consider the case of a vertical wall as shown in Fig 4. The vertical velocities are simply reversed across the wall. Since D = 0 in the fluid cell, it follows that the vanishing of D' is accomplished only if u' = + u1, in contrast to the requirement for u' for a free-slip wall.
In summary: (a) for a free-slip wall normal velocity reverses while tangential velocity remains the same; (b) for a no-slip wall normal velocity remains the same, while tangential velocity reverses.
For case (a) the pressure condition has been derived and has a simple form. For case (b) the no-slip wall, if vertical
φ' = φ1 ± gx δx ± (2ν u1/δx) (11)
where the sign is + for fluid to the left of the wall and - for fluid to the right of the wall. If it is horizontal
φ' = φ1 ± gy δy ± (2ν v1/δy) (12)
where the sign is + for fluid below the wall and - for fluid above the wall.
Velocity boundary conditions at a free surface are based upon the requirement that D ≡ 0 for surface cells. The easiest case to discuss is that of a surface cell which faces vacuum on only one side. Velocities for the other three sides are calculated in the usual manner, and that of the fourth side follows uniquely from the vanishing of D. Surface cells with two sides towards the vacuum are treated somewhat differently. For them we require δu/δx and δv/δy to both vanish separately; that is, each vacuum-side velocity is set equal to the velocity on the side of the cell across from it. This, of course, also assures that D = 0. A cell with three sides facing vacuum is relatively rare. The vacuum side opposite the fluid side is made to carry the velocity of that side; the other two vacuum sides which oppose each other are calculated to follow freely the effects of the body force and do not otherwise change. A surface cell with four sides towards the vacuum is similarly treated so that this isolated drop follows a free-fall trajectory.
The pressure boundary condition at the free surface is derived from the requirement of vanishing normal stress component [12] or of the equating of it to the applied external pressure. There are two reasons, however, why this condition is difficult to apply accurately. First, the normal stress components can be calculated only if surface orientation is known, and this can be accurately sensed in a finite difference representation only with great difficulty. Second, the various velocity derivatives at the surface are chosen on the basis of vanishing D, in a manner which apparently gives fairly accurate results in test problems, but which does not necessarily give the proper viscous stress in the surface skin. There are fairly complicated ways in which both of these difficulties can be resolved. For the test problems reported here, however, we have simply neglected the effect of the viscous stress in the surface boundary condition and equated surface pressure directly to applied pressure. The procedure will give valid results if the coefficient of viscosity is sufficiently small, and we have used computational experimentation as a means of testing the approximation. A more careful utilization of the correct boundary condition awaits further development of methodology.
Computational experiments have indicated considerable numerical stability for this computing technique. Of particular significance is the fact tlhat viscosity is not needed to insure stability. This is in contrast to many of the Eulerian techniques for fluid dynamics calculations in which the inclusion of an artificial viscosity term is required both for the treatment of shocks in compressible flows and for the elimination of stagnation fluctuations.
We have not succeeded in deriving the full necesary and sufficient conditions for numerical stability. Some necessary conditions for the case in which δx=δy are 4νδt/δx2 < 1 and (gdδt2/δx2) < 1. In the latter, g is the magnitude of the body acceleration, while d is the maximum depth. The first condition is the usual one for diffusional stability while the second is the incompressible-flow analogy to the Courant condition [12].
The question of accuracy has likewise been explored principally by experimentation. The succeeding section describes one set of calculations designed to test accuracy, and it is shown that detailed comparisons with physical experiments exhibit excellent agreement.
Since the computing technique is meant to be accurate for a variety of problems lacking crucial comparisons, it has been useful to develop an internal-consistency accuracy check. The natural choice is verification of overall incompressibility. For this purpose the computing code keeps a running account of the total number of cells which contain many marker particles. Santalo [14], quoting a result of Hadwiger, shows that if the marker particles were infinitely dense then the mean number of cells containing any should be
in which P is the perimeter length and A is the area. defined by the curve enveloping the particles. (We have dropped a term of order unity to account for each end of our perimeter usually being restricted to lie on a coordinate axis.) The first term dominates, while the second is a positive correction accounting for the fact that every partially filled surfuce cell is counted though totally filled.
With a sparse distribution of marker particles, the formula for N needs modification. Let λ ≡ (δsx/δx) ≡ (δsy/δy) be the ratio of particle spacing to cell size, here assumed the same in both directions. To account for λ ≠ 0, it is necessary to multiply the second term of Eq (13) by f(λ), a function for which an expansion is assumed to hold when λ is small: f(λ) = a0+a1λ _.... Now we require f(0) = 1-λ and f(1) = 0, leading to the result f(λ) = 1 - λ when higher-order terms are neglected. This lead, to the proposed formula.
In calculations for which it was verified that Di,j was negligible in every cell, we found excellent agreement between the observed value of N and that predicted by Eq (14), thus verifying the use of Eq (14) as an accuracy check. These tests examined both λ = 0.25 and λ= 0.50.
It should be emphasized that such a check is only a necessary, but not a sufficient, condition for accuracy. It could indicate good results, for example, even when the finite-diffence mesh was too coarse for proper resolution of the problem. Careful comparisons with the results of actual physical experiments give the only crucial tests of accuracy.
The flow dynamics behind a broken dam is a useful first application of the Marker-Cell computing method for three reasons. First, the experimental data available for comparison [15]. enable the technique itself to be proof tested. Second. the problem and its variants are ideal for demonstrating the versatility of the computing method as well as for showing the large amount of information that a calculation can produce. Third, the calculations demonstrate many features of this type of flow for which neither experimental nor theoretical information has previously been obtained.
We have performed calculations for two variations of the basic broken dam problem. In the first one the entire dam is instantaneously removed allowing the water to flow into a dry stream bed. In the second the dam is only partly broken, and the problem is equivalent to that of an impulsive partial opening of a sluice gate.
In addition to comparisons with experimental results, we have compared with those of shallow-water theory [16]. While the former showed very close agreement, the latter were found to be in considerable disagreement with the computer results.
All calculations were performed on the IBM 7030 (Stretch) Computer. The results shown in Fig 5 and Fig 6 were processed from the computer by the Stromberg-Carlson SC-4020 Microfilm Recorder. They are presented exactly as obtained from the computer with no retouching or other alterations. A typical calculation required from 10 to 30 minutes of computer time.
Fig 5 shows a typical sequence of marker particle configurations for the broken dam. Comparison of the results with those from the experiments of Martin and Moyee [15] showed excellent agreement in all respects. The deviations from shallow-water theory predictions are quite apparent. In particular, the front of the wave is much retarded in its initial motion and lacks the long sharp tip predicted by analysis. This last, incidentally, can be traced to the failure of the hydrostatic equilibrium assumption in shallow water theory rather than to the effects of viscosity described by Whitham [17].
The second type of broken dam problem reported here is equivalent to the instantaneous partial opening of a sluce gate. Several marker particle configurations are shown in Fig 6. In this case the grid of computational cells is completely omitted from the plot. There actually were twenty cells in the unit height of the initial water configuration.
One feature that is particularly well shown in Fig 6 is the tendency of lines of marker particles to spread in one direction as they compress in the other. It can even happen that internal computational cells may empty briefly, and it is necessary that a special technique be be employed to ensure that such cells continue to be treated as full of fluid. Since the marker particles do not enter into the calculations at all except to show free surface positions, the distortions of their appearance will then have no effect upon the validity or smoothness of the actual fluid dynamics.
Another computational feature is also indicated by the particle plots of Fig 6. Where the water is shallow near the front of its motion, the surface shows some steps. These can be attributed to the computational cells and can be removed by finer computational-cell resolution. Their effect is hardly appreciable even in the case presented here.
Fig 7 illustrates another type of plot available from the computer. It shows lines of constant pressure for the same calculation as that of Fig 6. Again, there is a computational feature requiring explanation. The zero-pressure surface boundary condition is applied at the centers of the surface cells, leading therefore to the inward dip of the zero pressure isobar at the opening at the earliest time. Note the crowding of the isobars as time advances and the initial large acceleration decrease.
A new computing method has been discussed which has proved successful for detailed calculation of the time-dependent, viscous , incompressible flow of a fluid with a free surface. The technique is applicable to studying a wide variety of problems in two-dimensional Cartesian coordinates and could easily be extended to other coordinate systems or three space dimensions.
Examples of problems already investigated are:
Calculations are planned to study the formation of waves by a linear explosion over the surface and the breaking of waves on a sloping beach.
Modifications to a cylindrical coordinate system will allow for more realistic studies of such problems as the effect of a spherical explosion and the splash of a spherical drop. The frequency of large amplitude waves is an interesting subject for either coordinate system. These and other modifications have been discussed in a recent report [20], that also presents numerous examples of the calculations.
This work was peformed under the auspices of the United States Atomic Energy Commission.
[1] B Alder, S Fernbach, M Rotemberg, editors Methods in Computational Physics, Academic Press, 1964, Vol 3.
[2] Y Mintz, University of California Report AFCRL 690 (1961).
[3] J O Wilkes, PhD thesis, University of Michigan (1963).
[4] J D Hell, H Z Barakat, University of Michigan, College of Engineering, Technical Report 1(1964).
[6] J W Deardorff, J Atmospheric Sci, 21, 419 (1964).
[7] C E Pearson, Sperry Rand Research Report SRRC-RR-64-17 (1964).
[8] D Greenspan, P C Jain, R Manohar, B Noble, A Sakurai, University of Wisconsin, Mathematics Research Center, Technical Summary Report 482 (1964).
[9] J E Fromm, A Method for Computing Nonsteady Incompressible, Viscous Fluid Flows, Los Alamos Scientific Laboratory Report LA-2910 (1963).
[10] J E Fromm, F H Harlow, Numerical Solution of the Problem of Vortex Street DevelopmentPhys Fluids 6, 975 (1963).
[11] F H Harlow, J E Fromm, Phys Fluids 7, 1147 (1964).
[12] L D Landau, E M Lifshitz, Fluid Mechanics, Addison-Wesley (1959), p51.
[13] F H Harlow, Los Alamos Scientific Laboratory Report, LAMS-2452 (1960)
[14] L A Santalo, Introduction to Integral Geometry, Herman and Cie, Paris (1953).
[15] J C Martin, W J Moyee, Phil Trans Roy Soc A244, 342 (1952).
[16] J J Stoker, Water Waves, Interscience Publishers, (1957) Chap 10.
[17] G B Witham, Proc Roy Soc A227, 399 (1955).
[18] F H Harlow, J E Welch (to be published).
[19] F H Harlow, J P Shannon, J E Welch, Liquid Waves by Computer, Science 149 1092 (1965).
[20] J E Welch, F H Harlow, J P Shannon, B J Daly, Los Alamos Scientific Laboratory Report LA-3425 (1965).