In this paper we discuss numerical methods for calculating multidimensional, transient, free surface flows interacting with general curved boundaries. To effectively model a free surface, three problems must be resolved: the surface must be numerically defined, a prescription must be provided to advance it in time, and appropriate boundary conditions must be applied at the location of the surface. The schemes chosen to meet these conditions must be compatible.
Before discussing methods that model free surfaces, we first review the basic notions of Lagrangian and Eulerian finite difference representations. In each representation the fluid region is subdivided into a finite grid, or mesh, of computing cells. Associated with each cell are local values of the appropriate field variables. The Lagrangian viewpoint utilizes a mesh of grid points embedded in the fluid and moving with it.
Thus, a free boundary will always be coincident with a mesh boundary, which simplifies the application of boundary conditions. The Eulerian viewpoint treats the mesh as a fixed reference frame through which the fluid moves.
In this case, some additional means must be supplied to locate a free boundary within the mesh. Furthermore, since the boundary is not likely to be coincident with the mesh lines, application of boundary conditions requires some interpolative subcell resolution scheme. For low amplitude free surface configurations, various methods are available that work more or less equally well. For example, several Lagrangian techniques have been developed to calculate incompressible flows with a free surface; see Hirt, et al. (1970), Brennen and Whitney (1970), and Chan (1975). These techniques, however, are not applicable to highly contorting or shearing free surface flows in which the Lagrangian meshes would be distorted so severely as to produce serious questions of accuracy.
Eulerian schemes are more frequently used for calculating free surface flows because they do not necessarily have the low amplitude restriction. For example, the Marker-and-Cell (MAC) method developed by Harlow and Welch (1965) is an Eulerian technique that has been extensively used to calculate multi-dimensional free surface phenomena. The basic technique has been improved by various users. In the original formulation the fluid is represented by marker particles, which are moved by the local fluid velocity. The fluid surface is defined by the location of these particles. The free surface boundary conditions are applied at the center of cells containing the interface. No attempt is made to determine the exact location of the surface within a calculational cell. This free surface treatment was improved by Chan and Street (1970). For low amplitude surface deformations, they defined the free surface by its height above the horizontal axis, which eliminated the need for marker particles. However, they also developed an alternative scheme that used marker particles only on the free surface. More importantly, in both schemes, they applied free surface boundary conditions at the actual location of the free surface. This practice eliminated small scale surface irregularities present in the original MAC formulation. Viecelli (1969) further extended the Marker-and-Cell method to include impenetrable, curved boundaries. He treats these boundaries analogous to a free surface, but applies a pressure at the boundary in such a way that the velocity of the fluid normal to the boundary vanishes. Thus, the difference between the treatment of cells containing free surfaces and those containing the impenetrable boundaries is that the cell pressures are computed from a Neumann (specified derivative) rather than a Dirichlet (specified value) condition. This technique was later extended by Viecelli (1971) to calculate flexible curved walls. Nichols and Hirt (1971) adapted and extended Chan's surface marker particle scheme to calculate highly contorted, multi-valued, free surfaces and applied better approximations for the normal and tangential stress conditions.
The assimilation of all these improvements into one general technique could result in a useful but very complicated technique. Realizing the need for these capabilities to be generally available in an easily used form, we have developed an Eulerian code, SOLA-SURF, which is capable of handling free and curved rigid surfaces that are single valued functions of height. The original SOLA code is based on the Marker-and-Cell method, but is written in a simplified form to facilitate its use by persons with little or no experience in numerical fluid dynamics. Because of its simplicity it is easy to add special boundary conditions to calculate, for example, flow past fixed or moving prismatic structures. In the SOLA-SURF version, which has the free surface capability, these structures may be completely submerged or may penetrate the free surface. The SOLA and SOLASURF codes are available from the Argonne Code Center and the basic algorithm and free surface scheme are described in a report by Hirt, et al. (1975).
Numerical techniques for calculating free surface flows in three space dimensions pose the same problems as the two-dimensional techniques discussed above, except they have the added requirement of making efficient use of computer time and storage. Nichols and Hirt (1973) developed a numerical technique for three-dimensional free surface flows. The technique is based on the Chan and Street version of the Marker-and-Cell method in which the free surface is defined by its height at the center of each vertical column of cells in the mesh. Because the surface must remain single valued and because of the manner in which boundary conditions are applied, the maximum surface slope must be less than that of a diagonal on the vertical face of any cell. While this does limit the amplitude of surface distortions, flow over and around some structures is possible in this technique.
Nichols (1973) has also employed marker particles for defining three-dimensional free surfaces. In this case, logic was established that permitted surface boundary conditions to be applied with arbitrary surface deformations. Although more general, there are problems associated with the use of marker particles for surface definition. The essential difficulty is associated with the lack of a suitable scheme to add additional surface markers in regions of surface expansion, for example, near the corner of obstacles that penetrate the surface.
Because of the difficulties associated with previous schemes, two new techniques are being developed. In one, surfaces are defined by a volume-fraction function that specifies the fractional fluid filled volume per unit volume. Surfaces are moved by fluxing the volume fraction from cell to cell. This eliminates the explicit setting and moving of a discretely defined free surface. The direction of maximum gradient of the volume-fraction function defines a surface normal direction, and this together with the volume-fraction values can be used to reconstruct an average interface location in each surface cell.
In the second alternative, a variable density fluid is assumed. The fluid has density normalized to unity, while in the void there is a very small density (for example, 10-6), No special boundary conditions are used in this method. For efficiency, these latter techniques are being developed initially for only two space dimensions.
All of the free surface schemes discussed in the following sections of this paper are couched in one basic solution algorithm, which is a direct extension of the Marker-and-Cell method. This algorithm is briefly described in section II. Following this a more detailed description, including the advantages and disadvantages, is given of free surface computational schemes that make use of the surface height function (section III), surface marker particles (section V), and the volume fraction and variable density schemes (section VI).
Before making a critical comparison of these different techniques, however, we illustrate how even the simplest method can be used to compute reliable and useful results. In section IV, following a description of the surface height function scheme, a study is presented of the added mass and damping coefficients computed for rectangular cylinders undergoing forced oscillation. These results are compared with the experimental data of Vugts (1968). An interpretation is also given of some nonlinear effects.
The finite difference solution of the Navier-Stokes equations is achieved by a direct extension of the Marker-and-Cell method. A stationary network of rectangular cells is used to divide the calculational region into a finite number of elements with which the fluid variables are associated. The primary field variables are the velocity components and the pressure. Each of the velocity components is specified at the center of the cell face to which it is normal and the pressure is specified at the cell center.
In several of the free surface schemes it is expedient to flag each mesh cell to denote whether it is an empty cell containing no fluid; a surface cell, which contains fluid but is not adjacent to an empty cell; or an obstacle cell, which is impervious to fluid flow. All remaining cells are interior fluid cells. Additional flags can be used to designate the rigid boundary cells, including interior obstacles, as either free-slip, which require zero normal velocity and zero shear stress, or no-slip, which require zero normal and zero tangential velocity at the boundary. Input and continuative outflow boundaries are also possible.
The fluid motion is numerically determined by advancing the fluid configuration through a series of small time increments. During each time step the solution to the momentum equation is obtained in two phases.
First, the known velocities and pressures from the previous time step are used to determine the fluid velocities in each cell, with the initial conditions used for the first time step. This explicit calculation does not necessarily ensure incompressiblity; therefore, in the second phase the tentative velocity field is adjusted through changes in the pressure field. This was accomplished in the original Marker-and-Cell method by solving a Poisson equation derived from the constraint that the velocity divergence be zero in each cell at the end of the time step. The basic idea of the present algorithm is to adjust the pressure in each mesh cell to drive the velocity divergence in that cell to zero. The pressure and velocity distributions must be obtained by iteratively adjusting these velocities in each cell in the mesh (see, for example, Hirt, et al. 1975).
The surface height function is generally more efficient and logically simpler than surface treatments involving marker particles. For this reason, this was the first method used to define the free surface in a threedimensional technique; see Nichols and Hirt (1973).
The fluid surface is single valued and is initially defined by specifying the surface height above the bottom of the mesh at the center of each vertical column of cells. The change in the surface elevation is determined by the local fluid velocity, that is, by the vertical component of the fluid motion plus the horizontal convection of the surface elevation from adjacent cell columns.
The free surface boundary conditions require that the normal and tangential velocities immediately outside the surface be chosen to ensure a zero transfer of momentum through the surface. Velocities normal to the surface are set to satisfy the incompressibility condition in the surface cells where the free surface is located. The tangential velocities in the cells immediately outside the fluid are set equal to the adjacent interior velocities. This is consistent with zero shear stress at the surface for an inviscid fluid. The pressure in surface cells is determined by a linear interpolation or extrapolation between the pressure in the fluid cell immediately below the surface cell and a specified pressure at the surface.
A new three-dimensional version of this free surface technique, SOLA- 3D, has been developed recently by Stein and Hirt (1975). The solution technique and the capabilities of the code are similar to the one described above, except that it has a variable mesh and the code has been written in a simpler, more versatile, and efficient manner. The following examples, however, have been computed using the original three-dimensional code.
All examples have been run on a CDC-7600 computer. A fast core memory of 64,000 words is available and a large core memory makes an additional 400,000 words accessible, although the large core is not a direct access core and does slow the calculation speed somewhat. The grind time, that is, calculation time per cycle per cell, for typical calculations varies from 0.8 to l.6 ms/cycle/cell. We have chosen two examples to illustrate the three-dimensional free surface features of this method. In Fig. 1, we show the perspective view of the free surface configuration resulting from flow over a rectangular structure. (Rigid bodies may be created in the fluid region by designating any number or configuration of cells as obstacle cells.) A 3 × 3 × 3 cell structure is located in a 15 × 15 × 7 cell mesh. The leftmost boundary is a fluid input boundary. The opposite end of the mesh is an outflow boundary. All other rigid boundaries are freeslip walls, that is, walls where the fluid is subjected to a zero normal velocity and a zero shear stress. The fluid is initially 3.5 cells deep and the Froude number is 2.0, based on initial depth and input velocity. The plots shown are at a nondimensional time of 2.0, after 20 calculational cycles. Perspective views from two separate viewing positions are shown. The plot on the left shows the fluid flowing through the mesh from left to right. The rightmost plot views the flow from a centered, downstream position.
Figure 2 shows two views of the free surface configuration resulting from fluid flowing from left to right past a partially submerged blunt body with a length to width ratio of 2.0, a width to draft ratio of 1.5, and a draft Froude number (based on the body draft and the input velocity) of 2.0. No frictional drag forces are assumed in the calculation. The wall boundary conditions are the same as in the previous example. At the calculation time of 3.5 (35 cycles) that is shown, the flow has reached steady state.
These examples are typical of the kinds of calculations possible using three-dimensional techniques with the free surface defined by a surface height function. The free surface is easily defined initially and the surface movement is straightforward. The free surface boundary conditions are easy to apply, even when adjacent to interior obstacles. The computer storage requirement for the free surface definition is minimal, since only two height values for each vertical column of cells are required, one old and one new time value. This free surface scheme can be used to calculate meaningful results in a low amplitude free surface environment, as demonstrated with the two-dimensional SOLA-SURF calculations presented in the following section. The only major disadvantage of this scheme is that it is restricted to a single valued free surface and, thus, to relatively low amplitude surface distortions.
To illustrate the utility of finite difference techniques we present a study performed to obtain the hydrodynamic coefficients of a rectangular cylinder. Added mass and damping coefficients are determined by computing the dynamic pressure force on the cylinder during forced heave. These calculations are compared with the experimental data of Vugts (1968).
The two-dimensional SOLA code, SOLA-SURF, which defines the free surface by a surface height function, was used for these calculations. The mesh size for the calculations varied from 100 to 250 cells in the horizontal direction, depending on the period of the forced motion, and there were 22 cells in the vertical direction. The mesh cells were of uniform size, δx= 0.1 and δy=0.2. Five cells were used to resolve the half width of the cylinder, which was located at a plane of symmetry at the left mesh boundary. To model the harmonic motion of the cylinder, special boundary conditions were added to the SOLA code. At the rigid boundary of the cylinder bottom, the cell pressure is derived from the constraint that the normal fluid velocity be equal to that of the body. The SOLA-SURF boundary conditions restrict the surface slope to be less than that of a cell diagonal. We eliminated this restraint at the vertical side of the cylinder, however, by aligning the side with a cell boundary line and adding boundary conditions that set to zero all velocities on and normal to the side of the cylinder.
The vertical displacement of the rectangular cylinder from its mean position is given by
a = a sin ωt
The linearized equation of motion for the cylinder in forced heave is, see for example Wehausen (1971),
(m + μ22) α..2 + λ22 α*2 + c22 α2 = Z where m = mass of the cylinder μ22 = hydrodynamic or added mass in heave λ22 =damping coefficient against vertical motion c22 = hydrostatic restoring coefficient against vertical displacement Z = external driving force
and the raised dot notation indicates differentiation with time. The total hydrodynamic force, that is, total force minus static force, Fb, acting on the body in the vertical direction is
Fb(t) = - μ22α..(t) - λ22 α* (t)
Assuming this hydrodynamic force is a harmonic function in time and can be written as
Fb(t) = γ sin (ωt - β) the hydrodynamic coefficients are obtained by equating these two expressions for Fb μ22 = (γ cos β) / (ω2 a) and λ22 = (γ sin β) / (ω a)
The amplitude γ and phase shift β are obtained by comparing plots of the displacement and dynamic pressure force F, acting on the cylinder as functions of time and measuring the shift in phase.
The calculated coefficients are compared with experimental data in Fig. 3. The coefficients are normalized by ρA and √(B/2g), where ρ is the fluid density, A is the mean submerged area, B is the cylinder beam, and g is the acceleration of gravity. The calculations were for B/T=2.0, where T is the cylinder draft at its mean location. The amplitude of motion normalized by B was 0.025 and 0.05. Vugts' experiments were conducted in fluid depths, normalized by B, varying from 4.50 to 5.625. The normalized depth was 4.0 for the calculated results. The calculated data generally agrees very well with the experimental data as shown in Fig. 3.
The discrepancy in the calculated damping coefficient at the normalized frequency of 1.25 is probably due to an error in determining the phase shift β, since it is very small at this beam width and frequency. To check this, another calculation was made, which is not shown in Fig. 3, for B/T=8 and ω√(B/2g)= 1.25. In this case the phase shift is much larger, and the added mass coefficient agrees closely with the experimental data. The calculated damping coefficient in this case is slightly lower than the experimental data, but is higher than the linear theory.
There is a discrepancy between the experimental data and linear theory for the added mass coefficient at normalized frequencies below 0.5.
Vugts explains this as being caused by experimental inaccuracies. With a normalized fluid depth of 4.0, the calculated coefficient at a frequency of 0.25 is slightly higher than the experiments. (This is the middle point marked (2) in Fig. 3 at this frequency.) However, calculations with a fluid depth of 2.0, marked (1), and a fluid depth of 8.0, marked (3), show clearly that the finite depth is the cause of the disagreement with linear theory. The added mass coefficient for the shallower fluid depth is slightly less than the experimental data, but with the normalized depth of 8.0 the calculated coefficient agrees with the linear theory, which assumes an infinite depth.
These calculations were done using a relatively coarse convergence criterion for solving the implicit pressure equation. Calculation times required from 3 to 15 minutes of CDC-7600 time, depending on the period of oscillation. At the low frequency, that is, ω√(B/2g)= 0.25, we observed a low amplitude, short wavelength disturbance superimposed on the hydrodynamic force. The disturbance appeared only on the positive increasing side of the force versus time curve, and can be eliminated by tightening the convergence criterion to increase the number of iterations per time cycle from 12 to approximately 60. The average magnitude and period of the force plots were not significantly changed by this, however, indicating that the coarser criterion is sufficient.
We observed some nonlinear effects in these calculations. As the rectangular cylinder is moving upward, a wake eddy forms near the bottom edge. On the downward motion this eddy dissipates somewhat as it is pushed from beneath the cylinder, but reforms as the cylinder velocity decreases, as seen in the velocity vector plot in Fig. 4. The development of this secondary flow is expected to increase the damping coefficient somewhat. According to Fig. 3, the calculated results are, in fact, higher than either linear theory or the experimental data. The experimental data at higher frequencies show significantly larger damping coefficients than predicted by the irrotational theory, but significantly smaller than the calculated results. We believe the calculated results are qualitatively correct, but are too large because they over-emphasize the secondary flow effects through numerical diffusion errors.
Another kind of nonlinear effect is responsible for asymmetries observed in the force versus time plot in Fig. 5. This calculation was at a normalized frequency of 1.25, with an amplitude of motion normalized by B of 0.095. At this larger amplitude, the damping coefficient is larger than the value calculated with an amplitude of 0.05, but the added mass coefficient is about the same. It is seen in Fig. 5 that the magnitude of the hydrodynamic force experienced by the cylinder at its maximum elevation is less than the force magnitude it experiences at its minimum elevation. Because the added mass coefficient is principally a reflection of the force amplitude, we can interpret this result as predicting a reduced added mass when the cylinder is at its crest and an enhanced added mass when at its trough. This result is reasonable as there is less fluid surrounding the body when it is raised, hence less that needs to be accelerated, and more when the body is at its low point. Linear theory does not predict this asymmetry because it applies free surface boundary conditions at the initial undisturbed boundary location.
We now consider a scheme to remove the low amplitude restriction in surface deformations. The original Marker-and-Cell method could be directly extended into three dimensions. However, the use of marker particles throughout the fluid volume would be costly and, as already mentioned, the application of boundary conditions at surface cell centers is sufficiently crude that bothersome fluctuations are generated. With these limitations in mind, Nichols (1973) replaced the surface height function technique described in section III with a free surface defined by marker particles. That is, particles are only used on the surface and boundary conditions are applied by interpolation at the surface position. The fluid motion is governed by conservative equations numerically solved with the solution algorithm described in section II.
Marker particles defining a free surface can initially be placed in any configuration within the three-dimensional mesh and are set in a staggered pattern for better surface coverage. For calculations involving inflow and outflow boundaries, particles are input as necessary to maintain coverage of new surface area flowing into the calculational region and are deleted as they flow from the calculational region.
A volume weighting scheme analogous to the area weighting scheme described by Welch, et al. (1966) is used to obtain interpolated cell velocities to move the surface marker particles. This is straightforward except in the vicinity of obstacles. For each of the three components of velocity, the volume weighting scheme must interpolate between eight neighboring velocities. Cells contiguous to these eight velocities are tested to determine if any are obstacle cells. If such is the case, the component of particle velocity parallel to the obstacle boundary must be determined as if the obstacle cell did not exist. If the component of particle velocity is normal to the obstacle boundary, the velocities used in the interpolation must not allow the particle to move into the obstacle cell.
In regions of large velocity gradients or expanding surface areas particles must be added to ensure uniform surface coverage. To add surface marker particles, an approximate surface plane in each surface cell is determined by first locating the approximate normal to the surface by considering the empty-full cell configuration about the surface cell. This surface plane is then projected onto the horizontal or vertical plane of the cell whose normal vector is closest in angle to the normal vector of the surface plane. This projected surface plane is then divided into 4 equal parts by the rectangular coordinate axes lying in that plane. One particle near the surface cell boundary is chosen from each of the three quadrants containing the most particles. The location of these particles is used to define the approximate plane of the free surface. Particles are added in this plane if the total number of particles in the cell becomes less than some prescribed number, such as half of the number used initially. This scheme functions very well, except in the regions downstream of obstacles and in flow around sharp corners. Without surface marker particles flowing into the cells in these areas, it is difficult to keep enough particles in the cells to accurately define the surface.
The free surface boundary conditions are imposed by setting the normal velocities to satisfy the incompressibility condition in the surface cells and the tangential velocities immediately outside the fluid to ensure zero shear stress. These are the same conditions used in the low amplitude technique described above and are appropriate for inviscid flow.
In the general case, however, there are 64 possible empty-full cell configurations about a surface cell, compared to 16 possibilities in the single valued surface case.
To proof test this technique, we made a variety of calculations, each of which tested a feature of the method. We have chosen two of these to present here. The grind time for the calculations varied from 2 to 8 ms/cycle/cell. Adequate free surface definition required from l6 to 36 particles per cell, which is 3600 to 8100 particles for a mesh with 15 × 15 cells in a horizontal plane.
The surface configuration resulting from flow over and around a submerged obstacle is shown in Fig. 6. The calculational mesh is composed of 15 × 15 × 8 cells and the obstacle size is 3 × 3 × 3 cells with the initial fluid depth equivalent to 4 cells. A fluid input boundary is at the left and an output boundary is at the right. All other mesh boundaries are rigid and free-slip. The Froude number, based on the input velocity and initial fluid depth, is 2.0. The surface plane in each cell is initially defined by 36 marker particles. The plots are at nondimensional times of 0.0, 0.5, and 1.0, and the calculational time step is 0.1. At each of these times, two perspective-view plots are shown, which are viewed from the same horizontal eye position but different vertical positions. In the plot on the right, the surface observed below the initial height shows the underside of the surface, while the surface above the initial height is the top side of the surface. The free surface and obstacles are transparent in these plots. The deformation of the free surface is caused by the initial impulse of flow past the obstacle and is completely formed by a time of 1.0. This surface deformation agrees qualitatively with the low amplitude calculation of a similar problem shown in Fig. 1.
We calculated flow past a partially submerged blunt body to test the free surface and particle movement boundary conditions near obstacles that penetrate the free surface. The calculational mesh is composed of 16 × 10 × 12 cells and the obstacle size is 10 × 2 × 8 cells located at a plane of symmetry. The draft is initially 2 cells. The fluid input boundary is at the left; the output boundary is at the right. All other boundaries are rigid, free-slip. Free surface configurations resulting from this flow are shown in Fig. 7. The plots are at nondimensional times 0.0, 1.5, 3.0, and 4.5, reading from left to right. The calculational time step is 0.1. The draft Froude number, based on the input velocity and obstacle draft, is 2.0. At a time of 1.5, the fluid has reached its maximum height on the obstacle front and has begun to fall back. At a time of 4.5, the fluid is beginning to fold onto the incoming fluid below.
In contrast to the surface height function, this technique is more complicated and less efficient, but does calculate multivalued free surface flows. The free surface is easily defined initially and the surface movement is straightforward, except near interior obstacles. More testing is involved in setting the free surface boundary conditions, because the free surface configuration can be more contorted than the low amplitude method. As described above, the major problem with the technique is in keeping the surface adequately defined in regions of highly contorted flow, especially on the downstream side of obstacle boundaries. This is a useful scheme for some multivalued free surface flows, but in its present form is not suitable in general.
Realizing the limitation of the surface height function method and finding the surface marker particle scheme not totally satisfactory, we attempted a method of defining the free surface through a volume-fraction function. This method has the simplicity needed in three-dimensional calculations and does handle highly contorted surfaces. This free surface scheme is used with the solution algorithm described in section II, except for the flagging scheme.
The cell centered volume-fraction function F is the ratio of fluid volume in each Eulerian cell to the cell volume, that is, F=1 in cells full of fluid, F=0 in empty cells, and F has intermediate values in surface cells. The surface is located at F=1/2. The basic equation solved is the convection equation written in conservative form for two space dimensions
∂F/∂t + ∂Fu/∂x + ∂Fy/∂y = 0
The F function is fluxed out of the positive face of each cell. For example, the F function is fluxed through the positive x-coordinate face of a cell by
Vx = δt u / δx
where δt is the time increment of the calculation, u is the velocity specified normal to the cell face, in the x-direction, and δx is the cell dimension along this axis. The sign of Vx defines the donor and acceptor cells (that is, the cells losing and gaining fluid volume, respectively). The amount of F fluxed out of this cell face each time step is
F = MIN{FAD |Vx| + MAX [(1.0 - FAD) |Vx| - (1.0 - FD), 0.0], FD}
where the subscripts A and D denote the acceptor and donor cells and AD is either A or D, depending upon the orientation of the interface, as explained below. The MIN feature negates fluxing more fluid than is available, while the MAX feature assures the fluxing of material from a partially filled cell. The amount of F fluxed is added to the acceptor cell and subtracted from the donor cell. Whether the acceptor or donor cell is used to determine the volume fluxed is dependent upon the orientation of the surface. The acceptor cell is used when the interface is convected normal to itself, and the donor cell is used when the interface is translating parallel with the interface orientation. However, if the acceptor is an empty cell, it is used to determine the flux, regardless of the orientation of the interface, thus allowing the donor cell to fill before any fluid is fluxed into this adjacent, downstream cell.
The values of the function F assigned to each cell serve as a flag, which is used in setting the free surface boundary conditions. The normal velocities are set to satisfy the incompressibility condition in surface cells. The tangential velocities immediately outside the fluid are specified to ensure zero shear stress at the surface. These are the same conditions set for the techniques described in section III and section V.
A useful calculation for testing this scheme, which was written as a two-dimensional code with plane or cylindrical coordinates and variable cell dimension capabilities, was one in which a cylindrical volume of fluid splashed into a quiescent pool. We defined the problem in cylindrical coordinates with the initial fluid volumes set by specifying values for the F function in appropriate cells. The cylindrical drop of fluid was composed of 7 vertical cells with a radius equivalent to 5 cells and the pool was 10 cells deep. There were 40×25 cells in the calculational mesh. A grind time of 2.6 ms/cycle/cell was required for this calculation. Many aspects of the code were tested with this problem. First, the calculation demonstrated its ability to calculate a body of fluid falling through space, while maintaining integrity of the body shape. Then, the free surface boundary condition logic and the volume-fraction function flux equations were given a severe test as the bottom surface of the cylindrical drop approached and collided with the free surface of the pool. Finally, momentum from the falling cylinder of fluid was imparted to the fluid in the pool as they coalesced, causing a contorted free surface, thus testing the codes' ability to handle these conditions. Figure 8 shows a sequence of free surface configuration plots for this calculation. Data is not available for direct comparison, but this calculation is in good qualitative agreement with a previous numerical calculation of a spherical drop of equal mass splashing into a shallow pool.
Because the free surface is defined by a cell quantity, this method is inherently less accurate than the other methods discussed. For very contorted free surfaces, setting the surface boundary conditions or selecting the correct F function to convect can become complex. This scheme is much simpler, however, than the surface marker particle treatment.
Attempting to further simplify the free surface treatment, we modified the equations of motion to treat a variable density fluid. The fluid density is normalized to unity, while in the void a very small density (for example, 10-6 ) is assumed. The donor-acceptor flux scheme is used to convect density in order to maintain sharp density interfaces. This scheme is presently being developed and does show definite promise. This scheme, if successful, will eliminate the explicit setting of the free surface boundary conditions, and will be useful for two fluid problems as well. A similar scheme has recently been developed by Spalding (1975).
The authors wish to thank R. K.-C. Chan for his interest and helpful comments on the material presented in this paper. This work was partially supported by the Office of Naval Research, ONR Task #NR 062-455 and by the United States Energy Research and Development Administration.
Brennen, C B, Whitney, A K, Unsteady, Free Surface Flows: Solutions Employing the Lagrangian Description of the Motion, Proceedings of the 8th Symposium on Naval Hydrodynamics, Pasadena, California, 117 (1970).
Chan, R K-C, Street, R L, SUMMAC-A Numerical Model for Water Waves, Stanford University Technical report No 135 (1970).
Chan, R K-C, Street, R L, , A Computer Study of Finite-Amplitude Water Waves, J Comp Phys, 6, 68 (1970).
Chan, R K-C, A Generalized Arbitrary Lagrangian-Eulerian Method for Incompressible Flows with Sharp Interfaces, J Comp Phys, l2, 311 (1975).
Harlow, F H, Welch, J E, Numerical Calculation of Time-Dependent Viscous Incompressible Flow of Fluid with Free Surface, Phys Fluids 8, 2182 (1970).
Hirt, C W, Cook, J L, Butler, T D, A Lagrangian Method for Calculating the Dynamics of an Incompressible Fluid with Free Surface, J Comp Phys, 5, 103 (1970).
Hirt, C W, and Cook, J L, Calculating Three-Dimensional Flows around Structures and over Rough Terrain, J Comp Phys, 10, 324 (1972).
Hirt, C W, Nichols, B D, Romero, N C, "SOLA - A Numerical Solution Algorithm for Transient Fluid Flows, Los Alamos Scientific Laboratory report LA-5852 (1975).
Nichols, B D, Hirt, C W, Improved Free Surface Boundary Conditions for Numerical Incompressible-Flow Calculations, J Comp Phys, 8, 434 (1971).
Nichols, B D, Hirt, C W, Calculating Three-Dimensional Free Surface Flows in the Vicinity of Submerged and Exposed Structures, J Comp Phys., 12, 234 (1973).
Nichols, B D, unpublished work performed for the Office of Naval Research under contract No. NA-onr-2-73, (1973).
Spalding, D B, personal communication (1975). 276
Stein, L R, Hirt, C W, SOLA-3D: A Solution Algorithm for Transient, Three-Dimensional Fluid Flows, Los Alamos Scientific Laboratory report in preparation (1975).
Viecelli, J A, A Method for Including Arbitrary External Boundaries in the MAC Incompressible Fluid Computing Technique, J Comp Phys, !t, 543 (1969).
Viecelli, J A, A Computing Method for Incompressible Flows Bounded by Moving Walls, J Comp Phys, 8, 119 (1971).
Wehausen, J V, The Motion of Floating Bodies, Annual Review of Fluid Mechanics, 3, 237 (1971).
Welch, J E, Harlow, F H, Shannon, J P, Daly, B J, The MAC Method, A Computing Technique for Solving Viscous, Incompressible, Transient Fluid-Flow Problems Involving Free Surfaces, Los Alamos Scientific Laboratory report LA-3425 (1966).
Vugts, J H, The Hydrodynamic Coefficients for Swaying, Heaving, and Rolling Cylinders in a Free Surface, International Shipbuilding Progress, 15, 251 (1968).