Particles are moved by applying forces and calculating the acceleration. Other attributes can change over time as well, which can enhance the representation of complex, dynamic objects.
Acceleration is then the rate of change of velocity over time:
Accelerations are calculated from forces applied to the particles.
The animation of a particle system using Newtonian dynamics can be summarized as follows:
for each time step calculate forces on all particles for each particle calculate acceleration:The variable f represents the force vector, and m is the mass. The integration shown above is known as Euler's method. It is the simplest and most common method, but is unstable, inaccurate and inefficient (Witkin and Baraff, 1993). Other higher order methods such as Runge-Kutta and midpoint are more accurate, and more stable. The system presented here uses the leap-frog method which will be discussed in detail in Section 3.3.integrate equations of motion:
Forces acting on particle systems commonly include gravity and damping:
Coupled particle systems are implemented by introducing binary forces which act on pairs of particles. A simple example is that of a spring:
A spring applies to a fixed connection between two particles. In contrast, forces can be defined which act over all particle pairs. In this situation particle interaction increases as two particles move closer together, and decrease as they move further apart. This idea stems from molecular dynamics, where long-range attraction forces and short-range repulsive forces determine the movement of molecules. This is commonly modeled using the Lennard-Jones potential energy function (Sprackling, 1985) which takes the form:
The first term models repulsion, the second attraction. Most literature reports using n = 12, m = 6, or smaller (Miller, 1989; Tonnesen, 1991). A and B vary according to the specific problem being solved.
The molecular dynamics approach has been useful in modeling fluid behavior (Miller, 1989; Terzopoulos, 1989; Tonnesen, 1991). A disadvantage to this type of force is that the calculation is O(N²), where N is the number of particles. This can be alleviated somewhat by limiting interactions to nearest neighbors, and through the use of creative data structures.
A collision is detected by checking the position of a particle against the plane made by the polygon in question. Witkin (1993) states:
Although general collision detection is hard, particle/plane collision detection is trivial. If P is a point on the plane, and N is a normal, pointing inside (i.e. on the legal side of the barrier,) then we need only test the sign of (X - P)*N to detect a collision of point X with the plane. A value greater than zero means it's inside, less than zero means it's outside (where it isn't allowed to be) and zero means it's in contact. (p. C11)A simple collision response involves modifying the normal component of the velocity vector. For a perfectly elastic collision (e.g. a golf ball) the normal component is simply negated. For softer collisions (e.g. jello) it can be multiplied by a constant r between 0 and 1. With r set to 0, the particle will not bounce at all. A contact force can be used to prevent a particle resting on the surface from being pushed into the surface:
which essentially cancels the normal component of the incoming force.
The idea of surface forces can be extended to the use of repulsive forces for preventing a particle from ever coming into contact with the obstacle. In essence a "force field" is set up around the obstacle by using a function which exponentially increases as the distance between the particle and the obstacle decreases (Terzopoulos et al., 1987; Monaghan, 1994).
Inserting a spring at the point of contact is another common collision response. This method is intuitive but computationally expensive. Stiffer springs require stiff equations and therefore require very small time steps (Moore and Wilhelms, 1988).
Individual particles can be rendered quickly and easily as point light sources or streaks of light. They can also be rendered using primitives such as spheres or cubes which can be drawn using either solid or wire frame models.
Surface rendering is a bit more complicated. Since a particle system does not specifically describe a surface, one must be computed. This can be done using implicit surfaces. Implicit surfaces are a class of curved surfaces which are defined as being the solution to some equation:
To generate a surface around a single point you can define a density field for that point. The field is generated using a function D whose value decreases exponentially as the distance from the point increases. An implicit surface would then be defined as the set of points where the density function is equal to some threshold T:
This method was introduced by Jim Blinn in 1982 who coined the term "blobbies" in reference to the resulting shape of the objects (Blinn, 1982). It has been used to model many types of deformable and soft objects (Wyvill et al., 1986a) including particle-based fluid (Miller, 1989; Tonnesen, 1991). The implicit surface can be rendered directly using ray-tracing or can be polygonized using well-known polygonization algorithms (Wyvill, 1986b; Lorensen and Cline, 1987; Bloomenthal, 1988).
The algorithm presented here renders the liquid model as a polygonal approximation of the surface defined by the particle system, with the option for individual particle rendering. The surface generation and polygonization are discussed in Section 3.5.
The mean pressure (P) on an area of fluid is defined as the ratio of the normal force acting on the area to the area. The pressure at a single point in the fluid is the limit that the mean pressure approaches as the area approaches zero (John, 1988). This can be demonstrated by looking at the force acting on a cubical element of fluid, as shown in Figure 1. Pressure can be exerted on an element by surrounding fluid elements or by the atmosphere if a surface is exposed to air.
The mean density (
)
of a volume of fluid is defined as the mass within a volume V divided
by that volume. The density at a point is the limit as the volume approaches
zero:
Density is related to the temperature of the fluid and the pressure exerted on the fluid elements: it will increase as pressure increases or as temperature decreases. The reaction to a change in pressure or temperature is different for a liquid than for a gas. When surrounding pressure is increased, a volume of gas will decrease substantially, causing a large change in density. When the same pressure change is applied to a similar volume of liquid, the volume will only decrease by a small amount, causing very little change in the density. The difference is due to the compressibility of the substance. A gas is much more compressible than a liquid, which means that its density can vary widely with pressure changes. The density fluctuations for a liquid are so small that a liquid is generally considered incompressible. The term "compressible fluid" will be used to refer to a gas, and "incompressible fluid" will refer to a liquid.
Fluid dynamics is the study of fluid motion in response to forces such as gravity and pressure. The equations of motion for a compressible fluid can be written in the form:
Equation 2.1 is referred to as the continuity equation. It describes the conservation of mass as the fluid evolves over time. Equation 2.2 is the momentum equation, which describes the conservation of momentum (mass multiplied by velocity), and equation 2.3 is the energy equation, describing the conservation of energy (John and Haberman, 1988; Sprackling, 1985).
In addition to these equations an equation of state must be used in order to fully describe the behavior of fluid. This equation defines a functional relationship between temperature, density and pressure (which varies according to the material being simulated). The form used for an ideal gas, sufficient to describe many compressible fluids, is:
where 
is the ratio of specific heats, a parameter which depends on the particular
gas being simulated.
A computational model of fluid involves representing the fluid as a set of small finite-sized elements, each representing a small volume of the fluid. The equations of motion can then be applied to the elements to simulate the behavior of the fluid over time. Tracking behavior and attributes as fluid flows past elements fixed in space is the Eulerian approach to describing fluid motion. Tracking the behavior of fluid elements as they move with the flow is the Lagrangian approach. Traditional CFD methods such as the Finite Element Method (FEM) and Finite Differences (FD) use the Eulerian approach. The work presented here is based on a Lagrangian CFD method known as Smoothed Particle Hydrodynamics (SPH). A brief comparison of the three methods, and a detailed description of SPH follows.
Smoothed Particle Hydrodynamics also represents fluid as a collection of elements with local fluid characteristics. But, instead of calculating flow from one element to the next, the actual movement of these elements in response to external forces is computed. Each one of these elements, termed a "particle", has a mass, position and velocity and is influenced by forces such as gravity, similar to the elements of the particle systems described earlier. However, an SPH particle also has local fluid characteristics such as density and temperature, and responds to pressure forces from surrounding particles.
The idea behind SPH is the determination of characteristics of fluid by interpolating from a set of non-ordered points: the particles. The interpolation is performed using a smoothing kernel W which is a weighted sum over particles within an area defined by a smoothing length h. As an example, the smoothed estimate of the density at particle i can be found using:
There are a number of possible choices for W , however the most advantageous (Monaghan, 1992) is:
Note that as q approaches 2 the value of W goes to 0 so that particles do not interact outside of two smoothing lengths. This fact improves the efficiency of the force calculations. Computing interactions among a restricted set of particles instead of across all particles can greatly reduce the number of computations per time step. Other methods of improving performance are discussed at the end of this section.
Computing the gradient of an interpolated value is done using the gradient of the smoothing kernel, giving a smoothed estimate of the gradient of the variable. The pressure gradient needed to compute the force on a particle is estimated using:
This is the manner in which the SPH versions of the equation of motion are derived:
The term 
is an artificial viscosity added to handle shocks. There are several variations
of this expression. The version implemented for this research is from (Monaghan,
1992):
Here 
is the mean density for particles i and j, 
and ß are viscous constants (usually set to 1 and 2 respectively), 
is a constant used to prevent numerical divergences, and c is the
mean sound speed.
The speed of sound for a particle i is:
which represents the speed at which sound travels through the fluid element represented by the particle.
Using these equations the particle positions and attributes are advanced using standard numerical integration methods such as leap-frog or predictor-corrector, and a variable time step. The size of the time step is controlled by the "Courant condition" which basically says that the time step should be smaller than the length of time it takes sound to travel over a specific distance. In SPH this distance is the smoothing length h. Section 3.3 describes the time step algorithm used in this research.
When updating the density, either equation 2.4 or 2.6 can be used. Most SPH calculations use equation 2.4 but equation 2.6 is more appropriate for incompressible flow (Monaghan, 1992; Monaghan 1994). The reasoning for this is discussed in section 2.2.3.
To put all of these equations in perspective, it may be helpful to summarize the steps a typical SPH simulation goes through. Particles are initialized in a configuration suitable to the problem being solved (e.g. evolution of matter orbiting the center of a galaxy). They are given an initial position, velocity, mass and energy and the system is evolved in a manner similar to the following pseudo code:
Update particle densities Update particle pressures while (time < end_of_time) for all particles Calculate acceleration due to pressure gradient (equation 2.7) Calculate rate of change of thermal energy (equation 2.8) for all particles Update position Update velocity Update thermal energy Update particle densities Update particle pressures Calculate new time step time += new_timestepSimulations in three-dimensions are usually run with 100 to 1000 particles. (Monaghan, 1992; 1994; Hernquist, 1989). The large number of particles combined with the inter-particle interaction algorithm results in a large amount of computations during every time step. This bottleneck can be alleviated by vectorizing or parallelizing the computations and/or using data structures which lend themselves to efficient nearest-neighbor traversal. Since particles only interact within two smoothing lengths the range of force calculations can be limited to "nearest neighbors." Algorithms based on linked lists or tree structures are commonly implemented for this type of computation (Theuns and Rathsack, 1993; Rhoades, 1992; Hernquist and Katz, 1989). A linked list implementation partitions particles into cubical cells which are one to two smoothing lengths wide. Interactions are then limited to particles in neighboring cells. Alternatively, particles can be stored in a hierarchical tree structure which encodes spatial relationships between particles.
In order to make a compressible fluid act as an incompressible fluid the pressure needs to be increased to a point where the density becomes so high that it is difficult to compress it any further. At this point small changes in density will result in very large changes of pressure, and conversely it will take a large amount of pressure to change the density by a small amount. The problem with such high pressures is that the Courant condition would now dictate extremely small time steps in the numerical integration.
Monaghan (Monaghan, 1994) uses this approach to adapt the SPH equations for incompressible flow. However, instead of approximating the real (incompressible) fluid with an exactly incompressible artificial fluid, he uses an artificial fluid which is more compressible than the real fluid. In other words, the simulated fluid is nearly incompressible, not entirely incompressible. This allows the time step size to be a little more reasonable, while still approximating a very stiff gas.
Incompressible SPH equations differ from traditional SPH equations in the density calculation, the velocity calculation and the equation of state. The standard density calculation (equation 2.4) is problematic because the density of liquid drops to zero at flow boundaries. If a smoothed estimate of that density is calculated using equation 2.4 the boundary particles will have a low density resulting in erroneous pressures from the equation of state. To alleviate this the particles are initialized with a given density and a rate of change for density is calculated based on the particles' relative velocity (equation 2.6)
The velocity of each particle is also adjusted to prevent particle interpenetrations and mixing. The XSPH variant (Monaghan, 1994; 1989) is used to coax particles to move with a velocity similar to the average velocity in its neighborhood which is said to be useful for high speed flow. It also helps keep particles of an incompressible flow orderly in the absence of viscosity. The variant is calculated using:
This variant is added to the velocity when advancing the particle position, but does not affect the velocity used for other SPH calculations.
Monaghan also uses a different equation of state, similar to that defined for water:
This equation produces large changes in pressure for small changes in
density, but if the system is in equilibrium (
= 0
) then the pressure becomes zero, which allows for larger time steps during
integration. (Note: During implementation of the SPH code presented here
equation 2.11 was found to be incorrect. See section 3.3)
Particle initialization and time stepping are the same as the original SPH (excepting the initial density mentioned earlier). When running a simulation Monaghan also uses damping to allow the system to settle into equilibrium before introducing perturbations (e.g. gravity or boundary movement) which initiate the motion to be studied.
