INESC, Rua Alves Redol, 9, 2º, 1000 LISBOA
In this article we present a method for the animation and visualization of turbulent fluid flow. The method is simple, fast and stable. It is based on well know methods from the field of Computational Fluid Dynamics, treating the fluid as a vorticity field. Vorticity is transported by a particle system. A uniform grid is also used to calculate a velocity field that moves particles to their new positions. Such a mixed method where free particles move over a fixed grid has very small computational costs, making it suitable for Computer Animation.
The method simulates the behavior of fluids in situations where the contact between fluid masses with different velocities gives rise to turbulence phenomena. It is suitable for the animation of gaseous fluids like smoke. Unlike previous algorithms, it is possible to generate turbulence over all scales, ranging from the macroscopic to the microscopic level. The algorithm is controlled by a small number of intuitive parameters, enabling animators to quickly take maximum advantage of it. The algorithm can also be parallelized easily owing to its particle nature.
The modelling and animation of fluids is one of the greatest challenges anyone working in Computer Graphics can face. Fluids are a constant presence in our lives, from the water that covers two thirds of our planet to the air that surrounds us. There is an understandable desire to portray such objects side by side with other more tractable geometric objects.
A polygonal mesh or a parametric surface are not adequate representations for a fluid because of its amorphous nature and irregular topology. A gaseous fluid tends to spread over all surrounding space, filling regions with variable density but with no definite border. Also common computer animation techniques provide realistic fluid motion only in the simplest of cases. Animations must resort to Fluid Dynamics laws to some extent if realistic results are desired. Computationally intensive algorithms are the penalty incurred by such approaches. Because of these difficulties in both the modelling and animation domains, fluids are one of the hardest real world objects to simulate on a computer.
In this article a method is presented for the animation and visualization of fluids exhibiting turbulent phenomena. Turbulence is a chaotic behavior of the fluid, characterized by fast variations of the fluid's velocity, both in space and time. Turbulence usually occurs when conditions of low viscosity and high speed gradients are present. A turbulent fluid can be visually identified by the presence of vortices. Such vortices continually form and evolve over time, giving rise to highly complex and appealing motions. Common examples can be a column of smoke rising through the air or the foam pattern left behind on the surface of water by a ship's propeller. From the point of view of both Computer and Traditional Animation, turbulence phenomena is frequently used to convey the notion of a fluid.
Several models for fluid turbulence have already been presented in the literature. In early works, turbulence was functionally defined as the superposition of band limited noises over consecutively smaller scales [Perlin85]. This ad-hoc definition of turbulence was used by [Ebert90] and [Saupe88] for the synthesis and animation of clouds. The results, although quite good for static images, were forcibly unrealistic when it came to simulating the dynamics of the clouds.
[Haumman88] was able to build arbitrary complex flow patterns by the addition of basic flow primitives. This primitives were simple solutions to the fluid mechanics equations for incompressible fluids and included uniform flows, sinks and sources of fluid and vortex effects. A static flow field could be built and particles would be advected through it.
[Sakas93] considered the turbulence of a fluid to be a stochastic noise superposed on a constant flow. The noise was synthesized from a power spectrum given by the Kolmogorov-Obukov distribution, which can be deduced from the equations of fluid flow [Landau63]. The turbulence patterns obtained by such spectral synthesis methods are everywhere uniform and isotropic. [Landau63] describes turbulence as a superposition of an infinite number of vortices, or eddies, with sizes varying over all scales. From the large scale eddies, energy is transmitted down to smaller ones without loss. The energy of the fluid is finally dissipated to the environment when it reaches the smallest eddies in the range of scales. Spectral synthesis with the Kolmogorov power spectrum takes this cascading of energies into account but fails to model the dynamics of the large scale eddies. Such eddies are responsible for the global, macroscopic, motion of the fluid and cannot be ignored if realistic results are expected.
[Stam93] also uses spectral synthesis to define a turbulent velocity field at the microscopic scale. To this, a deterministic flow field is added to describe the large scale motion of the fluid. A particle system is advected over the resulting field. The particles are treated as opaque blobs and rendered with a ray-tracing paradigm to give the fluid a gaseous look. The animations obtained with this method are fairly convincing but it is still not possible to model turbulence at large scales since the macroscopic flow field is deterministic. Also, energy and momentum transfers between this and the small scale velocity field are not considered. It is, therefore, not possible to model the onset of turbulence arising from the macroscopic motion of the fluid.
The model used in this work describes the fluid as a vorticity field. The vorticity is a measure of the amount of turbulence of the fluid at each point in space. If the vorticity field is known at some instant it is possible to obtain its evolution for all subsequent instants according to a given differential equation of motion. [Yaeger86] used this model to animate the atmosphere of the planet Jupiter. The vorticity field, corresponding to an initial cloud pattern, was sampled on a grid and the equation of motion was solved over that grid using pseudo-spectral techniques. Methods such as this, that sample and solve field equations over grids, are know as Eulerian methods. The grid must have a sufficiently fine resolution if the cloud pattern is to have a good level of detail. For a grid with a resolution of NN the method must solve a system of N2 differential equations, one for each node where a vorticity value is stored. Clearly, Eulerian methods are not particularly efficient when highly detailed animations are desired.
In opposition to Eulerian methods there is the possibility of describing the vorticity field as a particle system [Reeves83], this being the origin of Lagrangian methods [Leonard80]. A vortex particle (also know as a vorton in the Computational Physics literature) carries around a given amount of vorticity and can be thought of as an infinitely small whirlpool. When all these small whirlpools are added together the final turbulent pattern arises. In the method here described a grid is still used to keep an equivalent sampled version of the vorticity field (thus the method deviates from a "pure" Lagrangian formulation). A velocity field is calculated over the grid as the solution to a Poisson equation with suitable boundary conditions. The particles are then advected according to this velocity field.
To visualize the evolving fluid, particle positions are directly output to a scan-line renderer using a standard implementation of spot noise as described in [van Wijk91].
At present, the method is only applicable for two dimensional flow fields. It has, however, very low computational costs and can handle systems with large numbers of particles, which makes it an attractive method for the generation of complex animated textures. It also provides starting points for future extensions into the three dimensional domain.
The behavior of all kinds of fluids is governed by two equations: a continuity and a dynamics equation. For incompressible, inviscid, fluids and ignoring external forces like gravity these two equations take the form:
v(x,t) is a velocity field describing the fluid's flow for every instant, p(x,t) is the pressure distribution and is the density. Equation (1) represents a law of mass conservation. The amount of fluid entering any given volume must equal the amount of fluid going out. Equation (2) is the Fluid Mechanics equivalent to Newton's second law F=ma. The existence of a gradient of pressures inside the fluid gives rise to a flow v which seeks to eliminate it.
We will consider the flow to be purely rotational, i.e. there will always be a scalar field which allows the flow to be written as:
This transformation simplifies the dynamics since it can now be expressed with the field , known as the stream function, instead of a vector field. It can be seen by plugging (3) in (1) that a rotational flow automatically verifies the continuity equation.
At this point it is convenient to rewrite (2) as a function of vorticity in order to take turbulence effects into account. The vorticity (x,t) represents the amount of fluid circulation in every point in space and through time, being given by:
Taking that turbulence is really a superposition of vortices, one can see that the vorticity (4) is a good measure of the fluid's turbulence. Highly turbulent regions of the fluid will correspondingly have high vorticity values. Also note that in the case of a bidimensional fluid the vorticity is a vector always normal to the plane of the fluid. It therefore suffices to consider only its magnitude since the direction is fixed. By taking the curl on both sides of (2) one gets:
Notice that the curl applied on the left side of (2) managed to take the pressure gradient out of the picture. The dynamics are now further simplified since it only takes a stream function and a vorticity field to fully describe it. [Yaeger86] used this equation in his model to simulate the Jovian atmosphere. To try and find a more efficient method the variable is instead considered to be the vorticity carried by an infinitesimal fluid particle along its trajectory and through time. Equation (5) becomes:
The vorticity is therefore a conserved quantity of every particle. The equation of motion for a given particle i is now:
Simply stated, this equations means that each particle takes the velocity of the flow at its current position. Comparing (6) and (7) with (5) it can be seen that the dynamics of the flow are much simpler when written in Lagrangian coordinates than the equivalent version in Eulerian coordinates. One might already conclude that a particle system is a more adequate representation for a turbulent flow than a grid.
One equation is still necessary to have a complete vorticity model. This equation must relate the flow (in terms of its stream function ) to the vorticity field and is obtained by plugging (3) in (4) to arrive at the following Poisson equation:
It is important to note that such equation can only be obtained in the two dimensional case due to the orthogonality between and v(x,t). From the vorticity field at any given time it is now possible to get as a solution to (8) and then obtain the flow from (3).
A hybrid model is used where a particle system coexists with a fixed uniform grid [Hockney88], [Christiansen73]. The particles are responsible for the vorticity transport while the grid is used to solve the field equation (8). Each iteration of the algorithm can be briefly summarized in the following steps:
The algorithm is schematically represented in Fig. 1. This figure
shows the particle/grid duality and the circular flux of data
generated between the two domains. A particle system is advected
inside the Lagrangian domain while in the Eulerian domain bidimensional
fields are calculated, obeying differential field equations. The
transposition of any quantity from one domain to the other is
accomplished with a suitable interpolation scheme.
Figure 1. Hybrid algorithm.
The first step transfers the vorticity field from a Lagrangian
to an Eulerian representation. Figure 2 shows the vorticity interpolation
scheme. The particle vorticities are accumulated on the grid in
such a way that each particle contributes to the four grid points
that surround it. The particle vorticity Vp is distributed
through the grid points V1 to V4. A weighting factor
is attributed to each point according to the fractional area subtended
by the particle's position and the opposite grid point. For instance,
for V1 we have V1=A4Vp, for V2
V2=A3Vp, etc. In this way grid points nearer
the particle receive a higher percentage of its vorticity than
more distant points. This scheme also conserves the total vorticity
of the system.
Figure 2. Vorticity interpolation.
With the vorticity grid at hand the Poisson equation (8) is now solved relative to the stream function (x,y) on a rectangular domain a x b, c y d. Suitable boundary conditions must be imposed on the four edges delimiting the domain. Considering, for instance, the edges x = a and x = b one can have the following situations:
The laplacian operator in (8) is approximated with finite differences in order to solve the equation. The FISHPACK library was used, containing a set of Fortran routines to solve both Laplace and Poisson type equations over a variety of coordinate frames [Swarztrauber75]. Cartesian, polar and spherical surface frames are among the most used types of 2D frames.
At this point we have a discrete representation of the stream function over the grid. From this the fluid's velocity field v(x,t) can be obtained. For a given grid point (i,j) the curl (3) is approximated by the following differences:
where h is the grid spacing. For polar or spherical coordinate frames similar approximations to the curl operator are used. A particle's velocity is now calculated with an interpolation scheme similar to the one previously used to interpolate vorticities. Considering again Figure 2, where the grid points now store the velocities v1 to v4, the velocity vp of the particle is calculated through a bilinear interpolation:
During the advection stage, each particle is displaced according to:
where t represents the time increment. Equation (12) can be easily identified with the Euler method for the solution of ordinary differential equations. More robust methods, like the Runge-Kutta, could have been used, with increased computational costs.
The animation of the fluid is produced using a spot noise
approach to the rendering of the particles [van Wijk91]. The particle
system can be seen as a distribution of Dirac impulses in the
image plane. A Dirac impulse has an infinite frequency spectrum.
Trying to render it as is would inevitably lead to aliasing artifacts
in both the spatial and temporal domains.
Figure 3. Elliptic kernel.
To turn the particles into band limited signals a convolution is made with a gaussian kernel in the image plane. Considering a local coordinate space (x,y) for each particle where x is aligned with the particle's velocity vector (fig. 3), the kernel is given by:
where ws represents the spatial width of the kernel and wt is the temporal width. The rendering of the particle system is really a filtering operation with a low-pass, space and time variant, anisotropic filter (since it is not invariant with respect to plane rotations).
The convolution spreads the particle's influence over an elliptic area of dimensions x and y on the major axes. This spreading gives the image a "gaseous" look, being well suited for the rendering of gaseous fluids. A compromise must be met when choosing the filter's spatial width. Too small a value will not eliminate all high frequencies while a high value will lead to excessive blurring, causing an out of focus sensation on the viewer. As a good rule of thumb ws is set equal to the pixels width.
Temporal aliasing artifacts can be neglected when particles have small displacements from frame to frame. The parameter wt can be set to zero and the kernel becomes circular since the relation x = y now holds. The convolution with a circular filter can be performed more efficiently because there is no need to transform into particle coordinates as was previously done. The filter is rasterized one particle at a time in scan-line order, using incremental techniques and fixed point arithmetic. For a more detailed explanation refer to [van Wijk93].
The implementation of the algorithm has several speedup techniques. Grids have dimensions which are powers of two. In this way it is possible to index grid points with fast bitwise operators. Also the fact that particles are confined to move inside the grid allows us to use fixed point arithmetic to implement their advection. With 16 bit precision good quality animations could be produced with no significant errors.
Figs. 5a to 5c(12k) show frames from an animation of a turbulent jet. The complete animation is shown as an mpeg movie(35k). A lengthier animation can be seen here (127k). This kind of animation can be used to visualize the emission of smoke from a point source. Vortex particles are being created at a constant rate near the bottom of the images. They are uniformly distributed in the horizontal direction along the initial width of the jet. Particles are subject to a uniform flow field which drags them upwards at a constant speed. The turbulent flow field is calculated by the vorticity method over a 32x128 grid and is then superposed to obtain the final turbulence pattern.
Figure 4. Initial vorticity profile.
After some trial and error an initial vorticity distribution was chosen according to the profile in Figure 4, where 2W is the width of the jet at the base, along the horizontal direction. At first a uniform distribution in the range [-1,1] was tested but it turned out that particle vorticities cancelled each other in a random fashion and no significant turbulence phenomena ever developed. With the profile depicted in Figure 4 there are two opposing vorticity strips at each side of the jet. These strips are responsible for the large scale vortices that tend to form along the jet. The discontinuous transition in the profile is smoothed out after vorticity interpolation to the grid and forms an intermediate flow which gives origin to the arches connecting vortices on opposite sides of the jet. The vorticity of every particle is decreased at a constant rate in an attempt to simulate friction effects between the smoke and the surrounding medium.
Fig. 5d (4k) shows how motion blur can be achieved when the convolution filter (13) has a parameter wt 0. This can be used when the ejection speed of the jet is large compared to the speed of the shutter.
Figs. 6a to 6c (14k) depict the animation of a possible planetary atmosphere, much like the one modelled by Yaeger for the planet Jupiter. The atmosphere consists of several latitudinal cloud rings surrounding the planet. Vortex particles are constantly being introduced into the system of rings with a given initial vorticity. Each particle receives the color of the ring where it is initially placed. The concentration of particles per ring changes as a function of latitude. Equatorial rings show more complex dynamics and thus deserve a higher concentration of particles than polar ones where the flow pattern is rather dull. This change of concentration with latitude also accounts for the non-linearity of the latitude/longitude coordinate system.
The vorticity of each particle gradually decreases from frame to frame and the particle is removed from the system when its vorticity reaches zero. In this way the rings work very much like particle diffusers into the rest of the atmosphere. After some iterations the number of particles remains constant (the number of created particles equals the number of destroyed particles per iteration) and stable turbulence patterns generate between neighboring rings.
The model used here has several gross deficiencies from a scientific point of view, that became evident after the first simulations. The rings in the northern hemisphere circulate in a direction opposite to the rings in the southern hemisphere. Also, rings nearer the poles have higher speeds than equatorial rings. Several improvements include the modelling of the shear flow induced by the rotational motion of the planet and a more correct distribution of particle vorticities as a function of latitude.
The field equations were written in spherical coordinates and solved over a 64x64 grid. There is a periodic boundary condition on the longitudinal borders so the cloud pattern can wrap seamlessly around the spherical surface. Dirichlet conditions on the north and south pole borders force particles to circulate around without going over them.
The images are rendered in a two step process. A bidimensional color texture is first generated using the convolution technique for each of the RGB planes. Particles are convoluted with an elliptic kernel where the major axis (see Figure 3) is aligned with the latitudinal direction. The length of the axis is set to compensate changes in the concentration of particles. Rings near the poles have fewer particles and so need a good amount of latitudinal blurring to hide gaps in the resulting texture. In a second step, a standard scan-line algorithm renders the final images where the texture is mapped on a spherical object with a uniform surface color.
The vorticity method can not only be applied to problems with an intrinsically physical nature but also to a more broader class of Computer Graphics problems. Figs. 7a through 7d (50kb) show an example of image deformation based on an underlying vorticity field. The vorticity field is defined initially according to the image luminance. Bright areas of the image originate highly turbulent motion while darker areas tend to be passively deformed by the surrounding field. It can be seen that the bright area of the sky gives rise to a distinctly shaped large scale vortex. This direct mapping from luminance to vorticity is, of course, arbitrary. One could instead consider dark areas to be vorticity generators.
A vortex particle is created for every pixel in the image. Its color and vorticity are assigned according to the same attributes of the pixel. From here the vorticity algorithm starts working in the usual way. A swirling effect was achieved by including a flow field with the following expression in polar coordinates (r,theta):
The constants and define the twist and the rotational speed of a large vortical flow centered at the origin of the coordinate frame. The influence of (14) on the overall turbulent flow can be easily identified in the images.
To reduce both memory and computational costs the particles in the system were divided into two distinct groups of active and passive particles. A passive particle is advected much like everybody else but has no vorticity and consequently does not contribute to the system dynamics. It has a reduced attribute set, consisting only of its position vector and color. The vorticity field is defined only by the active particles which have the complete set of attributes (like vorticity and speed). The memory costs using this particle hierarchy are smaller for obvious reasons and also the time taken to interpolate vorticities to the grid is significantly reduced.
In the following table the times per frame are shown for each of the animations previously described. Ten algorithm iterations were performed for each frame in all cases. In the animations where the number of particles changes, values were shown for the frame with the larger amount of particles. The results were measured on a Silicon Graphics Indigo workstation with a R4000 processor running at 33 MHz.
|Jet||20 000||5,2 s||1,35 s|
|Atmosphere||192 132||37,2 s||47,6 s|
|Painting||(2) 327680||45,9 s||11,1 s|
The rendering time for the atmosphere sequence only measures the rendering of the bidimensional texture and not the final rendering of the texture mapped sphere. The discrepancy of this value relative to the rendering of the other sequences comes from the convolution with an elliptic kernel instead of a circular one. As already mentioned, a circular filter can be more efficiently scan converted.
It is easy to conclude the method has O(n) complexity, where n is the number of vortex particles. Particles do not interact with each other but only with the grid. The operations involved are advection, vorticity interpolation to and velocity interpolation from the grid and can be done in constant time regardless of particle distribution.
The time complexity for the Poisson solver is O(mn log(n)), according to the FISHPACK documentation, for a mn grid. It was verified, with the grid resolutions used in the animation sequences, that the computational cost of the Poisson solving operation was not significant relative to the cumulative cost of the other operations.
The algorithm can be easily parallelized. Particles are distributed to all the processors so that each one becomes responsible for the advection and interpolation of a small percentage of the total particle system. The grid is stored in a shared memory area and the field equation is solved with a conveniently parallelized version of the Poisson solver routine.
A method was presented for the simulation of turbulent behavior in fluids. The method is based on computational techniques, well established in the CFD domain for some decades. It is hoped this techniques may have taken some aspects of fluid dynamics out of the realm of large vectorizable supercomputers and have made them usable for medium sized workstations where they may find large application in Computer Graphics.
All the control variables of the model can be conveniently abstracted so that any animator, with no prior knowledge of Fluid Mechanics, can rapidly take advantage of it. For the animation of smoke, for instance, one only has to specify the concentration of particles, the thickness and initial velocity of the smoke and the desired degree of turbulence. A small set of tools can be implemented, each one aiming at a particular type of animation (smoke, clouds, turbulent textures) and with a reduced parameter set. Such computational tools could then be interfaced with commercial animation packages.
The algorithm is stable. [Hockney88] shows that the particle/grid method softens the short range interaction between particles. The vorticity of a particle is uniformly distributed inside the volume of a sphere centered at the particle's position. The radius of the sphere is a function of the grid spacing and the quantification errors made by the finite difference approximations. This fuzziness of the particle's vorticity makes the interaction force between particle pairs go smoothly to zero for short ranges. The possibility of hard collisions is avoided and the total energy of the system can be efficiently conserved.
The algorithm is also fast since its complexity changes linearly with the number of particles. Animation sequences using large numbers of particles were produced. These animations would have been quite impracticable or even impossible to produce with previous particle dynamics algorithms.
The main disadvantage of the method lies in its bidimensional nature. It can be used to produce animated textures which later can be mapped onto three dimensional objects. In a possible extension of the method to the three dimensional case, the vorticity field equation (8) would not take the simple form of a Poisson equation for which a variety of efficient and well tested routines exist. A special purpose numerical routine would have to be written and it is not known how this would compromise the stability and performance of the algorithm.