Reasoning about Fluid Motion I: Finding Structures

1 Introduction

It is commonly believed that there are two styles of sci­ entific thinking: analytical, a logical chain of symbolic reasoning from premises to conclusions, and visual, the holding of imagistic, analogue representations of a prob­ lem in one's mind so that perceptual and symbolic oper­ ations can be brought to bear to make deductions. Nei­ ther style is to be preferred a priori over the other. How­ ever, for problems whose complexity precludes a direct analytical approach, a certain amount of qualitative and visual imagination is needed to provide the necessary "feel" or "understanding" of the physical phenomena. Once the picture is clear, the analytical mathematics can take over and lead more efficiently to logical conclusions.

· The research was funded in part by NSF NYI Award ECS-9357773.

This "feel and physical understanding" is rather infor­ mal, imprecise, and apparently unteachable, but neces­ sary for scientists. My research goal is to formalize the visual style of thinking as computer programs, and to demonstrate the power of these programs by their abil­ ity to reason about the structure and motion of turbulent fluid flows.

The choice of fluid flows as a domain may seem ar­ bitrary and turbulence may add unnecessary difficulty to the project. However I believe the subject matter is fascinating, ripe for attack, more constrained than a commonsense theory of liquids, and, most importantly, is worthwhile because new tools for advancing turbulence research can have enormous scientific values.

Turbulent flows have been studied for many decades. But only recently has the capacity to perform direct nu­ merical simulation (DNS) of turbulent flows at moderate Reynolds numbers with enough accuracy been realized. 1 So for the first time detailed solution fields are avail­ able to scientists. Some of the important questions fac­ ing scientists are: How does one make observations from the data? How does one make theories based on new observations? How are data used to test theories?

Currently there is a large effort in Scientific Visual­ ization whose goal is to develop computer graphics to facilitate the presentation of large datasets. However the process of discovering interesting structures in the datasets and extracting physics from them is up to the human experts. Even with the help of modern visualiza­ tion software (such as AVS and Explorer), most human experts find the task time-consuming and prone to hu­ man visualization error. 2

1 The terminology is explained in sections 2 and 3.

In collaboration with experts in fluid dynamics from Yale and MIT, I am building computer tools that can reduce the post-processing time of CFD data by orders of magnitude. My approach is based on the following key idea: Active Visualization and Abstraction (AVA) . 3 By "active" I mean the reasoning process has three properties:

• Autonomous - the computer holds pictures in its mind so that perceptual operations can be used to make inferences which are otherwise too difficult to make by analytical methods alone.
• Purposeful- the visualization is done in the context of making observations and testing theories.
• Generative - the results of visualization guide the computer to make further observations.

By abstraction, I mean the process of reducing the en­ coding length of the datasets. A moderate-sized instan­ taneous flow field might consist of 0(10 7 ) grid points each with a 4-vector (three components for velocity and one for pressure) associated with it . It is hard to imag­ ine four numbers or even one at each point of space. Could we re-represent the flow field by semi-persistent objects that are much more compact and easy to visual­ ize? That the physical picture of moving and deforming objects, called vortices, is a useful level of abstraction for reasoning about fluids is an important claim of AVA.

What does AVA consist of? I conjecture that AVA consists of five core operations: 4

1. Aggregation
2. Classification
3. Re-description
4. Spatial inference
5. Configuration change

Aggregation is the grouping of primitive objects ac­ cording to some measure of similarity (e.g. , closeness , continuity, symmetry). For a familiar example, consider a vector field. Points in the vector field can be grouped into integral curves or orbits . Classification is the as­ signment of labels to the aggregate objects. Each la­ bel denotes a bundle of characteristic properties. Re­ description gives a more concise representation of the aggregate. These aggregates might possess characteristic shape and might touch or overlap. Deduction of these geometric relations is the province of spatial infer­ ence. The characteristic properties of an object might constrain what other objects can exist or change or de­ form in its neighborhood. Finally, the whole cycle of core operations can be reapplied at a more abstract level.

2 Further .discussion of this issue and current approaches in the visualization literature can be found in [Samtaney et al., 1994].

3 The concept is similar in spirit to the idea of "active vi­ sion" in the vision literature, which means the active control of camera movements and focus of attention to improve the robustness and stability of vision. I apply the concept to a higher-level cognitive process.

4 These core operations are distilled from a rethinking of how phase space analysis programs - like KAM [Yip, 1991] or MAPS [Zhao, 1994]- work. I suspect much of the work in mechanism analysis using configuration space (e.g. [Joskow­ icz and Sacks, 1991]) can be cast in the same framework.

A word about related works. Commonsense reason­ ing about fluids is a central problem in naive physics [Hayes, 1985b; 1985a]. The problem is hard because flu­ ids do not conveniently divide into discrete pieces that can be easily combined. Ken Forbus and his group have done important work in extending and partially imple­ menting Pat Hayes' ideas for representing fluids using both the contained-liquid and piece-of-stuff ontologies [Forbus, 1984; Collins and Forbus, 1987] . The work de­ scribed here is closer to the scientific end of the formal­ ization spectrum. As a consequence the theory is less general than Hayes' theory, but more relevant to the re­ search in the turbulence community. I expect the two lines of research to have fruitful interactions.

2 Properties of fluid motion

A fluid is either a gas or a liquid; it is not a solid. The behavior of fluid is interesting and often surprising. One can get a good understanding of fluids by thinking about the properties of water. Solids resist deformation ; fluids yield to any shear stress. Solids care how far they de­ form; fluids how fast they deform. The measure of the ease of deformation is viscosity. A more viscous fluid like honey deforms more slowly than a less viscous one like water. If the fluid is not moving, the only force it exerts is pressure, which always acts in the direction normal to any surface in the fluid . If the fluid is moving, a num­ ber of new forces come into play. Viscosity, an internal friction so to speak, causes a shear stress to develop be­ tween layers of fluid moving in different velocity. A fluid moving across a solid surface tends to drag the surface with it. This property is known as the no-slip condition, i.e., the velocity of the fluid is exactly zero at the solid surface. The frictional forces always act in the opposite direction of motion. A fluid moving over an asymmetri­ cal solid surface can create lift , a force perpendicular to the direction of motion.

Real fluids are quite complicated, so I will make some assumptions about the properties of the flow to simplify the discussion. The assumptions are convenient fictions, but in many interesting cases they are rather good ap­ proximations:

1. Fluids are incompressible. This is really the conser­ vation of mass. The volumetric flow rate is the same at every point in the fluid ; no fluid can accumulate anywhere.

2. Fluids are Newtonian. This gives a particularly simple linear relationship between the applied shear stress and the rate of deformation.

1. No temperature variation, Coriolis forces, nor elec­ tromagnetic forces exist that might affect the fluid. The Newtonian, incompressible fluid is already com­ plicated enough.

The payoff of these assumptions is large. The only relevant forces on the fluid, besides pressure, are inertial and viscous. The inertial forces keep the fluid going; the viscous try to stop it. The ratio between the two forces, a dimensionless number, is called the Reynolds number, denoted by the symbol Re.

One reason why the Reynolds number is an impor­ tant parameter is that the character of fluid motion is strongly dependent on it. A very high Reynolds num­ ber flow (Re ~ 1000) is dominated by inertial forces and tends to favor turbulence in which individual "fluid parti­ cles" move in a random unpredictable fashion even when the fluid as a whole is moving in a definite direction. Large viscous forces such as those in a flow with very low Reynolds number should damp turbulence and maintain a laminar flow in which fluid particles move more or less parallel to each other. In between the very high and very low Reynolds number flows are the transition flows whose character depends on the circumstances.

To get an idea of the different regimes of a real fluid flow, consider a steady incompressible flow past a sym­ metric "bluff body" (a non-streamlined body such as a circular cylinder or a sphere). See Fig. 1. At very low Reynolds number, when the flow is slow, the flow pattern consists of smooth symmetrical streamlines, the trajec­ tories of the fluid particles. As the Reynolds number is increased, a pair of symmetrical eddies or vortices ap­ pears in the rear of the cylinder where fluid particles curl around. If the flow speed continues to increase, the vortices become elongated and at some point break off, traveling downstream with the fluid. The phenomenon is known as flow separation. New vortices begin to form near the behind of the cylinder and shed alternatively into the cylinder wake. At sufficiently high Reynolds number, turbulence sets in and the wake begins to oscil­ late; the flow is unsteady and highly irregular. Vortices of many different length scales are also visible.

What really happens inside a turbulent flow? Is the motion of fluid particles purely random? Or is there organized motion within a background of smaller scale random fluctuation? Experimental evidence seems to fa­ vor the latter alternative. Nobody really knows for sure. But why do we care anyway? The reason is that although we do not understand turbulence, its effects are highly significant. Turbulence increases drag, mixing of materi­ als, and transport of heat. In designing an air transport, we might want to suppress turbulence in order to reduce drag; in a combustion engine on the other hand we might want to enhance turbulence to increase mixing rate. If large-scale structures do exist inside a turbulent flow and are found to be responsible for the enhanced mixing and transport properties, then it might be possible to control turbulence by direct interference with these structures.

3 Reasoning tasks

Given a numerical solution field of a turbulent flow, a spectrum of reasoning tasks can be defined. The follow­ ing list is roughly in the order of increasing complexity:

1. Making observations. Find out structures, if any, that exist in the solution field. Are there vortices? What are their shapes and sizes? How are they distributed? How are they created? How do they evolve and interact?
2. Making correlations. Determine how the shape and distribution of structures correlate with fluid ve­ locities, pressure, dissipation, and other statistical properties of the flow.
3. Incremental analysis. Given an instantaneous con­ figuration of structures, predict the possible short­ time behaviors.
4. Causal analysis. Explain and summarize the evolu­ tion of structures by a set of elementary interaction rules.
5. Testing theories. Given a hypothesis about struc­ ture formation or interaction, gather evidence to support or disprove the hypothesis.

4 Domain Theory

The entire theory of incompressible flow is contained in the N avier-Stokes equations:

The first equation expresses the conservation of mass; it is the incompressibility assumption. It says the veloc­ ity u , which has three components, has zero divergence. The name divergence is well-chosen for \7 ·u is a measure of how much u spreads out (diverges). If it is positive (negative) at a point P, then P is a source (sink) . Zero divergence means no source or sink inside the flow.

The second equation comes from Newton's second law of motion; it equates fluid acceleration with applied forces . The acceleration terms, the left hand side of the second equation, represent the change of velocity in time and space. On the right hand side, the first term -\7 p is the pressure gradient. The second term, vis the kine­ matic viscosity which is the ratio of the viscosity and density of the fluid; the product of v and the second spatial derivatives of the velocity \7 2 u is the force due to viscosity.

Since we can associate with every point (x, y, z) in space and instance t with a velocity u, we call u(x, y, z, t) a velocity field. A field is any physical quantity that takes on different values at different positions and time. u( x, y, z, t), like a magnetic field, is a vector field be­ cause at each space-time point , we associate a vector with three components. The pressure field p(x,y,z,t) like temperature is a scalar field - one number for each point.

While useful for some specialized flows, the velocity field and its topology are not a particularly good repre­ sentation of turbulent flows. The velocity field not only is not Galilean invariant (what an observer sees depends on how fast she is moving), but also can change quite unpredictably from one instance to another. As a con­ sequence, we cannot visualize what is happening - not easily. Just as some scientists prefer to think in terms of interaction of charges and magnets rather than the re­ sulting magnetic field, many fluid dynamicists prefer to think in terms of a vector quantity called vorticity, de­ noted by the symbol w, which is the curl of the velocity, w= \7 x u. The curl of u at a point P measures how much the vector field u curls around P.

Suppose you float a paddlewheel on a bathtub. If it starts to turn, then the point it is placed has a non-zero curl. A region with a large curl is an eddy, a whirlpool. The curl of u is a vector; its direction is assigned by the right hand rule: if the water surface is the xy-plane and the paddlewheel turns counterclockwise, the curl points to the upward z-direction.

The reason we introduce the vorticity field w(x, y, z, t) is that it simplifies the description of fluid motion and gives ontological primitives which are easier to reason with. A vortex line is an integral curve of the vorticity field. The set of vortex lines passing through a simple closed curve in space is said to form the boundary of a vortex tube. Vortex lines and vortex tubes have nice invariant properties. In particular, they can be treated as material objects. The truth of this last statement fol­ lows from the so-called Helmholtz laws of vortex motion. For incompressible, inviscid (i .e., v = 0) flow, Helmholtz proved the following theorems [Batchelor, 1967]:

1. The vorticity field has zero divergence (because the divergence of a curl is always zero).
2. Vortex lines move with the fluid, i.e., fluid particles that at any time lie on a vortex line continue to lie on it .
3. The strength r of a vortex tube, defined as circula­ tion r = fs w · ndS is the same for all cross-sections S of the vortex tube and is constant in time.

The first and second theorems explain why vortex lines or tubes can be treated as material bodies. The third theorem is just an expression of the conservation of an­ gular momentum. The skater's spin is a nice illustration of the third theorem at work. By bringing in her arms and thereby shrinking the cross-section, the skater spins faster because the total angular momentum is conserved.

Turbulent flows are not inviscid, but at very high Reynolds number the viscous effect is very small except in a thin region close to the :;olid boundaries. So it is reasonable to expect the Helmholtz laws to be approxi­ mately correct.

Although the exact analysis of the motion of a con­ figuration of vortex lines and vortex tubes can be very complicated, the simplest situations can be understood by a few qualitative principles:

Qualitative Rules of Vortex Motion

1. A vortex line accelerates velocity on one side and slows down on another. A lift force perpendicular to the vortex line is generated due to the Bernoulli effect.
2. A vortex line is convected by the fluid which exerts a drag force to push the vortex line in the direction of fluid motion.
3. A vortex tube stretched (compressed) in one direc­ tion increases (decreases) the velocity components in the other two directions.
4. A bent vortex line, conceptualized as a space curve, exerts a self-induced motion along its binormal and the effect is largest at points of maximum curvature [Arms and Hama, 1965].
5. Vorticity can only be created at velocity discontinu­ ities and solid boundaries.

1. When two antiparallel vortex lines are brought close together, they can break apart and reconnect. (See Fig. 2)
2. When a vortex line meets a free-slip boundary, it can reconnect.

The first four rules are consequences of Helmholtz laws; they are good approximations even when applied to vor­ tex tubes with finite core size as long as the core is rel­ atively thin. The fifth is a kinematic consequence of the no-slip condition and the conservation of vorticity. The last two rules describe how the topology of vortex lines can change. The detail of the reconnection mechanism is still an open issue, but reconnection appears to happen experimentally.

The significance of these qualitative rules is that they allow the synthesis of the velocity field - at least the rotational component of the velocity. Main qualitative features of the flow field can often be deduced without complicated numerics. The use of these rules in incre­ mental motion analysis will be the subject of a sequel paper.

5 Automatic Extraction of Vortex Structures

5.1 Aggregating vortex lines

A vortex line is the basic building block of a vortex struc­ ture. Because of the divergence-free property, a vortex line, like a magnetic field line, does not start or stop in the interior of the fluid; it tends to run in a closed loop. However, a real turbulent flow always has a background of randomly fluctuating vorticity. So it is reasonable to assume that only the relatively high-intensity vor­ tex structures remain coherent from one time instance to another. To distinguish a "coherent" vortex from a mathematical vortex, I propose the following definition for aggregation:

Definition A coherent vortex is a compact bundle of adjacent, high-intensity vortex lines that are geometrically similar.

There is some degree of arbitrariness in any definition of coherent structure because the notion of coherence is informal. I believe the definition here is more faithful to the mathematical definition of a vortex tube. A coherent vortex can be tube-like or sheet-like depending on the shape of its cross-section. Unlike a mathematical vortex, a coherent vortex can start or end in the interior of the fluid.

Aggregating vortex lines to form coherent structures is not straightforward. Previous researchers [Moin and Kim, 1985; Robinson, 1991] have found that vortex lines are sensitive to initial conditions. Nearby vortex lines can diverge rapidly. If the initial conditions are not cho­ sen carefully, the resulting vortex lines are likely to re­ semble badly tangled spaghetti wandering over the whole flow field, making the identification of organized struc­ ture extremely difficult . This might explain why vortex lines have not been widely used for structure identifica­ tion. I believe the search algorithm below is the first successful structure identification based on vortex lines.

Given a numerical vorticity field, the search algorithm finds all coherent vortices. The key idea in the algorithm is the adaptive control of the cross-section of the vor­ tex tube: the cross-section is shrunk (expanded) when the vortex lines on the boundary of the cross-section are converging (diverging). The algorithm has the following steps:

1. Find all grid points that are local extrema of vortic­ ity magnitude and greater than a threshold. These are the seed points.
2. On the plane normal to the largest vorticity vec­ tor component at the seed point, find an isocontour centered at the point. The contour, discretized into points, represents the initial cross section of the sur­ face.
3. Use the advancing front method (explained later) to interleave the advancement of the cross section by integration and the tiling of the surface.
4. Use the geometry of the tiles to decide shrinking or expanding the cross section locally.
5. The wavefront is periodically adjusted globally by computing its convex hull, and the diameter and width of the hull.
6. The forward integration terminates when the circu­ lation on the cross section falls below certain thresh­ old.
7. Reconstruct the surface by integrating the last cross section backwards until it reaches the initial cross

section. No wavefront adjustment is needed in this step.

8. Remove the weak vortex lines.

The advancing front method is implemented as fol­ lows. (See [Hultquist, 1992] for details.) The vortex surface consists of a list of ribbons. Each ribbon has two tracers: a left and a right . As its tracers are advanced, the ribbon is tiled by triangular meshes in such a way to keep wavefront nearly perpendicular to the integration direction. The aspect ratio of the quadrilateral formed by the last pair of alternating left and right triangles in each ribbon is used to control the local adjustment of the wavefront.

Vortex lines can twist and turn, so the cross sections can get distorted quite a bit as the surface is developed. The vortex-finding algorithm keeps track of the number of tracers. If the number exceeds a threshold, which is indicative of many highly distorted quadrilaterals, the entire cross section is replaced by the convex hull of the projected wavefront on a plane normal to the vorticity vector at the centroid of the cross section. The convex hull is useful for three purposes: (1) it reduces the num­ ber of tracers, (2) it re-orders the tracers into adjacent positions along the vertices of an oriented polygon, and (3) important shape information of the cross section such as its diameter and width can be computed in linear time from its convex hull [O'Rouke, 1994].

5.2 Aggregation results

As the test case, I use the DNS results of a free sur­ face turbulence provided by Professor Dick Yue in the Ocean Engineering Department of MIT. The turbulence is generated by a shear flow in a 128 3 rectangular box with periodic boundary conditions in the x and y direc­ tion. A 4th order polynomial interpolation is used to compute the interpolated vorticity vector from the grid values, and an adaptive 4th order Runge-Kutta integra­ tor to integrate vortex lines. Hundreds of structures have been constructed by the algorithm. A typical result is shown in Fig. 3a. The vortex is reconstructed from 10 vortex lines. The computation including rendering (done by AVS 5) takes about 90 seconds real time on a Spare 10/51.

The algorithm is not sensitive to the initial choice of isocontour value: it is self-adjusting. Contrast this with an ordinary integrator. The vorticity lines obtained di­ verge and are tangled (Fig. 3b). Moreover, small changes in the initial isocontour can result in drastically different vortex line patterns.

5.3 Classification and Re-description

The reconstructed vortex must be interpreted in order to perform spatial inferences and incremental motion anal­ ysis. Classification is the assignment of labels to vortices. The assignment is determined by the shape of the vor­ tex. A vortex can be tube-like or sheet-like. Tube-like vortices are further differentiated by the curvature func­ tion along their axes. It is important to identify the local curvature extrema because these extrema cause a self-induced motion (section 4).

To obtain a concise representation of a bundle of vor­ tex lines, a vortex is re-described as a generalized cylin­ der (GC). A GC consists of a spline, a cross-section, and a sweeping rule [Binford, 1990]. The GC representation of a vortex is computed by the following steps:

1. Pick the vortex line with the highest integrated vorticity as the candidate spline.
2. Compute a scale-space representation of the cur­ vature of the candidate spline [Witkin, 1983].
3. The stable local extrema are chosen as knot points.
4. At each knot point P, compute where the vortex lines intersect the plane at P with the plane nor­ mal equal to the vorticity vector at P.
5. Fit an ellipse to the intersection points to obtain a cross-section of the vortex.
6. The spline of the GC is obtained by spline-fitting the centers of the elliptical cross-sections.
7. Compute the intrinsic shape descriptions of the spline, i.e., its curvature and torsion.

By varying the criteria for stable curvature extrema, one can obtain generalized cylinders of different resolu­ tions. Fig. 4 shows a rather fine generalized cylinder representation of the reconstructed vortex first shown in Fig. 3a.

6 Conclusion

This paper presents three novel ideas:

1. A zeroth-order formalization of the visual style of thinking as a cycle of five core operations: aggrega­ tion, classification, re-description, spatial inference, and configuration change.
2. A formalization of a theory of fluid flow based on new ontological primitives: vortex line, vortex tube, and coherent vortex, and a list of qualitative inter­ action rules.
3. A new vortex-finding algorithm based on vortex lines.

The implementation of spatial inference (such as de­ termining spatial relations among vortices) and config­ uration change (such as incremental analysis of vortex motion) will be the subject of a sequel paper.