|
[In this technical article originally printed in Game Developer magazine, Neversoft co-founder Mick West looks at how to efficiently implement fluid effects - from smoke to water and beyond - in video games, with example code.]
Fluid effects, such as rising smoke and turbulent water flow, are everywhere in nature but are seldom implemented convincingly in computer games. The simulation of fluids (which covers both liquids and gases) is computationally very expensive.
It's also mentally expensive, with even introductory papers on the subject relying on the reader to have math skills at least at the undergraduate calculus level. In this two-part article, I will attempt to address both these problems from the perspective of a game programmer who's not necessarily conversant with vector calculus. I'll explain how certain fluid effects work without using advanced equations and without too much new terminology.
I'll also describe one way of implementing the simulation of fluids in an efficient manner without the expensive iterative diffusion and projection steps found in other implementations. A working demonstration in ZIP form with source code accompanies this article; example output from this can be seen in Figure 1.
FIGURE 1 Smoke output is shown from the accompanying code.
Grids and Particles
There are several ways of simulating the motion of fluids, but they all generally divide into two common styles: grid methods and particle methods. In a grid method, the fluid is represented by dividing up the space a fluid might occupy into individual cells and storing how much of the fluid is in each cell.
In a particle method, the fluid is modeled as a large number of particles that move around and react to collisions with the environment, interacting with nearby particles. Let's focus first on simulating fluids with grids. The simplest way to discuss the grid method is in respect to a regular two-dimensional grid, although the techniques apply equally well in three dimensions. At the most basic level, to simulate fluid in the space covered by a grid you need two grids: one to store the density of liquid or gas at each point and another to store the velocity of the fluid. Figure 2 shows a representation of this, with each point having a velocity vector and containing a density value (not shown). The actual implementation of these grids in C/C++ is most efficiently done as one-dimensional arrays. The amount of fluid in each cell is represented as a float.
The velocity grid (also referred to as a velocity field, or vector field) could be represented as an array of 2D vectors, but for coding simplicity it's best represented as two separate arrays of floats, one for X and one for Y.
FIGURE 2 The fluid density moves over a field of velocities, with a density stored at each point.
In addition to these two grids, we can have any number of other matching grids that store various attributes. Again, each will be stored as a matching array of floats, which can store factors such as the temperature of the fluid at each point or the color of the fluid (whereby you can mix multiple fluids together).
You can also store more esoteric quantities such as humidity, for example, if you were simulating steam or cloud formation.
|
In graphic application one can use a faster (second order in time) trapezoidal rule for the integration :
x_i - x_i* = (u(x_i, t+dt) + u(x_i*, t)) / 2
This is an implicit equation for x_i* (the position back in time) which could be solved by estimating u using the 2 first term of the Taylor expansion around x_i
x_i* = x_i - dt ( 1 + dt/2 Grad u(x_i,t)^-1 * ((u(x_i,t) + u(x_i, t+dt)/2)
where all terms on the RHS are computed at grid point x_i.
The biggest problem with Stam algoritm is the use of linear interpolation, which is not appropriated for rendering of fluid. One early drawback of semi-Lagrangian schemes was that, since the basis functions of common high-order interpolants such as cubic Lagrange or Hermite are not positive definite, a scalar that is initially bounded between 0 and 1 will not necessarily remain so. However, this problem can be easily fixed by a local clipping, as proposed by Bermejo
and Staniforth (1992). This local clipping produces a quasi-monotone scheme. The term quasimonotone means that the scalar field can develop wiggles, but it cannot exceed its range at the previous time step. This is different from traditional clipping, where under-shootsand over-shoots are reset to minimum and maximum background values chosen, based on physical considerations, at the beginning of the computation. As shown by Bermejo (2001) this clipping is equivalent to a linear combination of a higher order interpolant and a low-order (monotone) Lagrange interpolant. To guarantee (discrete) conservation of scalar, two approaches are proposed by Priestly (1993) and Bermejo and Conde (2002). Personnally I prefer the one from Priestly (1993), since it can be used with the quasi-monotone interpolation scheme previously described. The complete reference is
Priestly, A. 1993 A Quasi-Conservative version of the Semi-Lagrangian Advection Scheme. Mon. Wea. Rev., 121
http://www.google.com/patents?id=ySgIAAAAEBAJ&printsec=abstract&zoom=4
see http://ams.allenpress.com/perlserv/?request=get-abstract&doi=10.1175%2F1520-0493
(1996)124%3C2046%3AMFFSLT%3E2.0.CO%3B2
Using this stability condition for the choice of the timestep, any two trajectories cannot cross each other. In this case, there is no need to keep a list of points of origin and mass fraction.
ass%20project).zip
from my understanding this code is light enough to perform well on mobile devices, I was happy to read the possible improvements suggested in this article! :)
Du/Dt =du/dt + (u.grad)u= + nu grad^2 u
it dosent matter in the end as the "vorticity confinement, and projection" steps are effectively taking care of it, but the physical reason they are needed is the pressure!
Du/Dt=du/dt + (u.grad)u= ## -1/rho grad P## + nu grad^2 u
However, there are a few problems which arise later - particularly dealing with the decompression step. A special case exists wherein a thick wave gradually gains speed - the cells before it cause the wave to move faster towards them - since the wave is more than one cell thick, it isn't balanced off by the rear. As it moves to a faster and faster velocity, it begins to move across multiple cells at the same time - resulting in an ever increasing velocity which is not slowed in the slightest.
This causes multiple glitches in the system - as soon as a thick wave starts, it will grow and fill the entire fluid (as the decompression algorithm will continue to spread it out in a given direction.)
My quick fix - dissalow cells that have more than a given amount of fluid to decompress - the backward advection will pull them out, and they will eventually slow down enough to spread out to a normal state.