Free Surface Fluid Flow
Table of Contents
Fluid flow problems often involve free surfaces in complex geometry and in many cases are highly transient. Examples in hydraulics are flows over spillways, in rivers, around bridge pilings, flood overflows, flows in sluices, locks, and a host of other structures. A capability to computationally model these types of flows is attractive if such computations can be done accurately and with reasonable computational resources. To be useful, simulations should be much faster and less expensive than using physical models.
Many computer programs can solve the partial differential equations describing the dynamics of fluids. Not many programs are capable of including free surfaces in their simulations. The difficulty is a classical mathematical one often referred to as the free-boundary problem. A free boundary poses the difficulty that on the one hand the solution region changes when its surface moves, and on the other hand, the motion of the surface is in turn determined by the solution. Changes in the solution region include not only changes in size and shape, but in some cases, may also include the coalescence and break up of regions (i.e., the loss and gain of free surfaces).
In this note a computational modeling technique for fluid flows with arbitrary free surfaces is discussed. The technique is based on the Volume-of-Fluid (VOF) technique. This technique has many unique properties that make it especially applicable to flows having free surfaces. The goal of this discussion is to show why the VOF approach offers a natural way to capture free surfaces and their evolution with great efficiency.
A good recommendation for the VOF method is to demonstrate its capabilities on a simple hydraulic flow problem, one that is far from trivial. The example selected is of flow over a step. This flow has conceptual simplicity and good experimental data available for validation (see N. Rajaratnam and M.R. Chamani, “Energy Loss at Drops,” J. Hydraulic Res. Vol. 33, p.373, 1995).
Prototype Hydraulic Flow with Free Surfaces
Figure 1a shows the flow problem after it has reached a steady-state condition. The overflow (sheet of liquid or nappe) leaving the top of the step has both an upper and lower free surface. At the bottom of the overflow a pool has formed between the overflow and the face of the step, while downstream, liquid is flowing to the right with a flat, steady surface. Strictly speaking, the flow conditions in the pool region are not steady because turbulent mixing is generated in the pool by the impinging fluid. There is, however, an average configuration and that is what is reported in the experiments.
For all practical purposes the flow is two-dimensional, that is, it does not have any significant variation in the direction normal to the illustration in Fig. 1a. In actuality, to have an air space above the pool there must be some opening to the atmosphere otherwise it would close up.
The flow speed at the top of the step is critical, that is, it has a speed equal to or greater than the speed of surface waves, so that no disturbances from downstream can penetrate through this region to affect flow upstream (to the left of the step), which is why the flow is exceptionally smooth and steady in that region.
There are many geometric features in this problem that can be compared with a numerical simulation; such as flow heights before and after the step, the angle of the overflow stream when it strikes the bottom and the depth of the pool formed under the overflow. Additionally, an important comparison for practical applications is the amount of energy (i.e., kinetic plus potential) lost by the flow in passing over the step.
Simulation of Prototype Problem
Figure 1a is from a simulation. For this example all of the geometric and material properties used in the experiments were used in the simulation. The height of the step used in the laboratory test is 62cm and the fluid is ordinary water (density=1.0 gm/cc and dynamic viscosity=0.01dynes/cm). The depth of water entering the computational region was 15.5cm and was given a near critical velocity of 123.0cm/s. Of course, gravity was in the vertical direction with magnitude g=-980cm/s2.
Because some turbulence was expected to develop in the pool to the left of the overflow, a turbulence model (the Renormalization Group or RNG model) was used in the simulation. Subsequent simulations without a turbulence model produced very similar results, which is not too surprising since most of the important elements of the flow are smooth (i.e., non-turbulent) inflow, overflow and outflow streams.
The simulation region shown in Fig. 1b is 170cm wide and 100cm high and has been subdivided into a grid of equal sized rectangular cells consisting of 80 cells in the horizontal direction and 60 cells in the vertical direction, for a total of 4800 cells. This grid is used as the basis for finite-difference approximations of the governing differential equations of fluid dynamics (the Navier-Stokes equations). The number and size of the grid cells was chosen with the goal of capturing the smallest expected features of the flow. The number can be easily increased or decreased if the results seem to warrant some adjustment. In fact, it is often a good idea to repeat a simulation with a change of resolution to make sure that the solution is not too sensitive to such changes.
The left boundary was a specified velocity boundary (also with a specified fluid height). The right boundary was an outflow boundary where all flow quantities have a zero gradient normal to the boundary to encourage a uniform outflow. The top and bottom boundaries are rigid walls, while in the third direction the boundaries were treated as planes of symmetry (i.e., walls with zero viscous drag). The surface of the step was also treated as a free-slip boundary.
Initial conditions could have been set to roughly approximate the expected flow arrangement, but since the flow configuration is one of the things that one would like to compute, especially for situations where one doesn’t know what the distribution of fluid is likely to be, a simpler approach is needed. Because a transient flow simulator was used for this example a simple initial condition could be defined that consisted of just a block of fluid on top of the step, Fig. 1a with the same horizontal velocity and height assigned to the left boundary. The simulation then followed the development of the steady flow, which occurs after about 8.0s. The simulation was run out to a time of 10.0s to assure that steady conditions had been reached. Figure 2 shows two intermediate times; 2.b at 0.2s and 2.c at 0.5s plus the final time in 2.d at 10.0s.
It should be noted that what starts as a single, connected free surface changes to two independent free surfaces (upper and lower nappe surfaces) after the fluid strikes the bottom. No difficulties are experienced with this separation of the flow into portions flowing to the left and right of the impact point on the bottom boundary. This will be discussed at further length in the next section.
Comparisons between experiment and simulation are given in the following table and are in excellent agreement.
| Experimental Results | Simulation Results | |
|---|---|---|
| Outflow Height/Step Height | 0.094 | 0.094 |
| Pool Height/Step Height | 0.41 | 0.41 |
| Angle of Nappe at Bottom | 57° | 59° |
| Energy Loss/Initial energy | 0.29 | 0.296 |
In view of these results it might be expected that a considerable amount of computational time would be required to achieve such accuracy. In fact, the total cpu time on a desktop Pentium 4, 3.20GHz computer was only 88s. Such a short computational time requires explanation and that is the purpose of the following sections.
Why the VOF Technique Works Well
There are a few general concepts about computational methods and the VOF technique in particular that can be used to gain an understanding of how and why VOF works so efficiently.
Basic Theory
All numerical methods must use some simplification to reduce a fluid flow problem to a finite set of numerical values that can then be manipulated using elementary arithmetical operations. A typical procedure for approximating a continuous fluid by a discrete set of numerical values is to subdivide the space occupied by the fluid into a grid consisting of a set of small, often rectangular “bricks.” Within each element an averaging process is applied to obtain representative element values for the fluid’s pressure, density, velocity and temperature.
Simple equations can be devised to approximate how each element’s values interact with neighboring elements over time. For instance, the density of an element can only change when there is a net flow of mass exchanged between an element and its neighbors (i.e., conservation of mass). The material velocity that carries mass between elements is computed from the conservation of momentum principal, usually expressed in the form of the Navier-Stokes equations, which uses the pressures and viscous stresses acting between neighboring elements to approximate the changing fluid velocities in the elements.
This idea of an element interacting with its neighbors is essentially what is meant by a partial differential equation; that is, evaluating the effects of small changes caused by the variation in quantities nearby. Partial differential equations are typically derived in engineering text books as the limit of approximations made with small control volumes whose sizes are then reduced to infinitesimal values. In a numerical simulation the same thing is done except that the control volume sizes cannot be taken to the limit because that would require too many elements to keep track of. In practice, the goal is to use enough elements to resolve the phenomena of interest, and no more, so that computing times are kept to a minimum.
Arithmetical operations associated with an element generally involve only simple addition, subtraction, multiplication and division. For instance, the change of mass in an element involves the addition and subtraction of mass entering and leaving through the faces of the element over a fixed interval of time. A simulation requires that these operations be done for thousands or even millions of elements as well as repeated for many small time intervals. Computers are ideal for performing these types of repetitive operations very rapidly.
Simulating fluid motion with free surfaces introduces the complexity of having to deal with solution regions whose shapes are changing. A convenient way to deal with this is to use the Volume of Fluid (VOF) technique described next.
The VOF Concept
The VOF technique is based on the idea of recording in each grid cell the fractional portion of the cell volume that is occupied by liquid. Typically the fractional volume is represented by the quantity F. Because it is a fractional volume, F must have a value between 0.0 and 1.0.
In interior regions of liquid the value of F would be 1.0, while outside of the liquid, in regions of gas (air for example), the value of F is zero. The location of a free surface is where F changes from 0.0 to 1.0. Thus, any element having an F value lying between 0.0 and 1.0 must contain a surface.
It is important to emphasize that the VOF technique does not directly define a free surface, but rather defines the location of bulk fluid. It is for this reason that fluid regions can coalesce or break up without causing computational difficulties. Free surfaces are simply a consequence of where the fluid volume fraction passes from 1.0 to 0.0. This is a very desirable feature that makes the VOF technique applicable to just about any kind of free surface problem.
Another important feature of the VOF technique is that it records the location of fluid by assigning a single numerical value (F) to each grid element. This is completely consistent with the recording of all other fluid properties in an element such as pressure and velocity components by their average values.
The VOF Technique
For accuracy purposes it is desirable to have a way to locate a free surface within an element. Considering the F values in neighboring elements can easily do this. For example, imagine a one-dimensional column of elements in which a portion of the column is filled with liquid, Fig. 3. The liquid surface is in an element in the central region of the column, which will be referred to as the surface element. Because we assume the values of F must be either 0.0 or 1.0, except in the surface element, we can use this to locate the exact position of the surface. First a test is made to see if the surface is a top or bottom surface. If the element above the surface element is empty of liquid, the surface must be a top surface. It the element above is full of liquid then, of course, the surface is a bottom surface. For a top surface we compute its exact location as lying above the bottom edge of the surface element by a distance equal to F times the vertical size of the element. A bottom surface is similarly located a distance equal to F times the vertical size of the element below the top edge of the surface element. Locating the surface within an element in this way follows from the definition of F as a fractional volume of liquid in the element.
Calculating surface locations in one-dimensional columns is simple, accurate and requires very little arithmetic. In two and three dimensional situations, however, computing a location is a little more complicated because there is a continuous range of surface orientations possible within a surface cell. Nevertheless, dealing with this is not difficult. A two-dimensional example, Fig. 4, will illustrate a simple way to not only compute the location of the surface, but also to get a good idea of its slope and curvature.
As in the one-dimensional case, it is first necessary to find the approximate orientation of the surface by testing the neighboring elements. In Fig. 4 the outward normal would be closest to the upward direction because the difference in neighboring values in that direction is larger than in any other direction. Next, local heights of the surface are computed in element columns that lie in the approximate normal direction. For the two-dimensional case in Fig. 4 these heights are indicated by arrows. Finally, the height in the column containing the surface element gives the location of the surface in that element, while the other two heights can be used to compute the local surface slope and surface curvature.
In three-dimensions the same procedure is used although column heights must be evaluated for nine columns around the surface element. Although a little more computation is needed, it consists primarily of simple summations in the columns and then sums and differences of column heights for evaluating the slope and curvature. Based on this discussion, the reader should now see how the fractional fluid volume can be used to quickly and easily evaluate all the information needed to define free surfaces.
There are two remaining issues to deal with. One issue is that a simulation like that in Figs. 1 and 2 is only solving for the fluid dynamics in regions where there is fluid. This is another reason for the computational efficiency of the VOF method. The region occupied by fluid in the flow over a step problem is much less than half of the open region in the computational grid. If it were necessary to also solve for the flow of gas surrounding the liquid, then considerably more computational time would be required. In order to perform solutions only in the liquid, however, it is necessary to specify boundary conditions at free surfaces. These conditions are the vanishing of the tangential stress and application of a normal pressure at the surface that equals the pressure of the gas.
A second issue is that movement and deformation of a free surface must be computed by solving for the fraction of fluid variable, F, as it moves with the fluid. Because the variable F is discontinuous (i.e., primarily 0.0 or 1.0) some care must be taken to maintain this discontinuity as it moves through a computational grid. In the VOF method, special advection algorithms are used for this purpose.
Illustration of Free-Surface Tracking by VOF Technique
Figure 5a is an illustration of how well this works; the fluid volume fraction is colored uniformly in each grid element to represent its value in that element. The free surface is sharply defined nearly everywhere. Only in the lowest and narrowest part of the nappe is there any noticeable loss of a sharp fluid fraction distribution, Fig. 5b. This was expected because in this region the nappe is less than three elements in thickness and this allows some of the smaller F values associated with partially filled surface elements to mix in with the central element, which should have a value of 1.0. For computational purposes this doesn’t really matter because the simulation method treats elements interior to the liquid as though they are pure liquid elements.
It should also be pointed out that in the region shown in Fig. 5b turbulence and air entrainment are observed in actual experiments. Thus, the appearance of fluid fraction values a little less than unity is somewhat realistic. This is not entirely accidental because the intersection of jet of liquid with a pool, which is responsible for turbulence and air entrainment, is also responsible for the “entrainment” of fluid fraction values into the interior of the liquid.
Summary
At first it may seem somewhat magical that a computer can simply perform repeated arithmetic operations on arrays of numbers and produce a realistic simulation of a complex, time-dependent, fluid dynamics problem. It was the purpose of this discussion to explain an approach that does this with relatively elementary procedures.
Using a simple, but non-trivial, hydraulic flow example it has been demonstrated that computational simulations can produce detailed results in excellent agreement with physical measurements. It has been further demonstrated that the simulation, which was based on the Volume of Fluid (VOF) technique, uses simple approximation methods that are both accurate and efficient.
Clearly, real world examples involving complex hydraulic structures such as those used in hydroelectric power stations, must consume more than the few seconds of computational time used in our example to obtain useful results. Nevertheless, those results can be generated in reasonable times (both man and computer) and contain a richness of detail rarely possible in physical experiments. For examples visit our water and environmental application pages. In addition, the ability to easily test the influence of just about any kind of change in geometry, flow condition or fluid property is another powerful reason to employ simulations. Current software and hardware for hydraulic flow simulations offer a significant cost advantage over traditional physical modeling.
Postscript
The first detailed description of the VOF method was in 1981 by C.W. Hirt and B.D. Nichols, J. Comp. Phys., 39, p.201. All simulations appearing in this article were performed with the commercial software package FLOW-3D developed by Flow Science, Inc. This program uses an enhanced variant of the VOF concept called TruVOF.
