VIRTUELLE WELTEN

A PARTICLE SYSTEM WITH CONSTRAINTS

Paper by Arno Schödl, Virtual Worlds II project, TU Berlin, summer term 1997

A particle system is a cloud of particles, each of these has, among other arbitrary attributes, a certain mass, a location and a velocity. They move according to forces acting on them, either external forces like gravity or drag, or forces depending on other particles, like spring forces or gravity between planets.

Particle systems are not only good models for obvious physical phenomena like dust, snow or water, but they also work for chains, oscillating springs or rigid and not-so-rigid bodies, as long as we restrict the particles´ movement in a suitable way.

Given a force acting on a particle we update the particle´s position x according to:

,

where v is the particle´s velocity vector, p is its position vector and m its mass. Note that I omit any vector designator. The vector F is the force acting on the particle. It is important to notice that F can change with time or any quality of another particle, so it´s only accurate at a definite point of time within a particular environment, that means positions and attributes of particles.

When we store all the n particles´ positions and velocities living in 3D space in one vector x, we have a 6n dimensional vector. We can now form an ordinary differential equation, describing the behavior of this system state vector:

f determines for each component of the vector its derivative, for position components it just copies the according velocity component, for velocity components it calculates the current force in this direction, and returns the acceleration.

Assume we somehow got f, that is the derivative for each component of our state vector x. How do we numerically solve the differential equation?

NUMERICALLY SOLVING DIFFERENTIAL EQUATIONS

Imagine a ball swimming on a lake with some complex currents, with the ball always moving with the current. Now imagine the lake has not two, but 6 times the number of particles dimensions (I know, you can´t!), then you have the picture of the state of our particle system moving in phase space along the vector field of derivatives.

At the beginning, the ball is at , the initial state of our particles. Our task is to take samples of some reasonable points of the lake, determine the current there by evaluating f and then decide where the ball will be carried.

The first technique is very obvious, it´s called Euler´s Method:

It assumes the current in the lake does not change at all, the ball is carried along with the current stream sampled at the position of the ball. This method is very prone to the illustrated accuracy and stability problems:

In the taylor expansion of the state vector x at the time

all but the first two terms are omitted. But, there are ways to use higher terms of the taylor expansion. The problem is that our function f provides only , but no higher derivatives. By evaluating f at more points though we can approximate the higher terms. The greater number of evaluations required to do a step are more than compensated for in terms of computational costs by the longer steps possible with these methods.

The update rule using one term more is called midpoint method, its update rule is as follows:

a) Do a normal euler step:

b) Evaluate f at the midpoint of the step:

c) Take the step using the midpoint value

Even using higher derivatives is possible, the general method is called Runge-Kutta. But, although we can reduce our computational costs using these methods, we want to visualize the movements of our particles, so it´s no use taking too large steps because we want to get a smooth animation. A compromise has to be made between frame rate and the accuracy and stability of our particle system differential equation evaluation.

FORCES ACTING ON PARTICLES

The next step is how to define our function f giving the derivatives of velocity and position. The derivatives of the positions are the velocities, a copy operation is all that is to do. The derivatives of the velocities depend on the forces acting on the particles. Here are some useful force functions:

GRAVITY

g is the gravity acceleration, a vector pointing down with the magnitude of the gravity (on earth 9,81 m/s²).

VISCOUS DRAG

k_d

is the drag coefficient, viscous drag describes the force acting on a particle moving in a viscous liquid like honey.

SPRING FORCE

where F_a and F_b are the forces on the respective particles, l = a – b, r is the rest length of the spring, is the derivative with respect to time, k_sis a spring constant, and k_d is a damping constant.

What does this formula do? The value in the square brackets calculates the value of the spring force, the direction is given by the vector –l= b – a, going from a to b, moving particle a closer to b. The value of the force consists of two components, first there is the spring force, is attracting force is proportional to difference between the current length of the spring and its rest length. Second there is a damping force, resisting movement of the spring like viscous drag and therefore dampening the spring´s oscillations.

CONSTRAINTS


Particle system get a lot more powerful if we can define conditions for legal positions of certain particles, fixing particles somewhere in space, fixing the distance between particles, allowing them only to move on a predefined path and so on. One is tempted to implement these constraints as very firm springs. But remember what we said earlier about the differential equation evaluation. It works well for smooth and slow changes in. Firm springs oscillate very fast, therefore requiring very small time steps to be simulated correctly, consuming a lot of computation power, or instabilities will occur, destroying our whole simulation.

The problem is that spring forces caused by illegal positions result in accelerations, these in turn result in velocities one time step later, and another time step later positions are changed with these velocities to enforce the contraints, so using the Euler method it takes at least two time steps to correct a violation. Our goal is to simply inhibit forces that lead to constraint violations, enforcing the constraints at each time step without delay.

How do we generally define constraints? We define m functions c_i, one for each of m contraints, which have to be zero when and only when the respective constraint is met. The function must be differentiable. The reason for this will become clear later.

Our functions need the position of all particles as an input, so we define a state vector q that this time only has the positions as components, not the velocities, so in 3D for n particles it has 3n components.

We do the same thing with the velocities, getting the derivative of our position with respect to time, that is . We define the vector-to-vector-function C simply as the vector of all the outputs of the constraint functions:

We assume that the initial positions are legal with all contraints, so . The velocities are legal, too, at least for the moment, so Cdoes not change with respect to time, that means . In order to keep these two conditions valid from one time step to the next we have to keep , which we will do by introducing new forces that eliminate the forces which violate constraints. Assume we do nothing against constraint violation, we just calculate our force vector, let´s call it f, as we did before. What will happen to C, or more specifically, to ?

To describe changes of C when parameters of C are changed we introduce the Jacobian matrix J of C:

What does this mean? The functions c_i get a vector as input, so they represent higher dimensional landscapes, with the height of the landscape being the value of c_i. In its m rows the Jacobian matrix J contains the gradients of the constraint functions c_i, i=1...m, having as many columns as there are components in the input vector q.

Now we derive C(q) twice with respect to time, using the chain rule:

, where , that is the matrix of gradients of C derived with respect to time, and , the acceleration calculated from the masses and the accumulated forces of the particle.

This equation describes the change of C we have to counteract. The first term describes the change of C because of the particles´ velocities. At the moment all positions are legal, the particles are in a c_i = 0 area of the landscape described by the contraint functions c_i and all velocities are legal, that means, all the landscapes are flat at the point where we are. But maybe ahead of us we start to move uphill? That is what the first term checks for and applies enough force to counteract. I think it is obvious now that the constraint function must not have a step, or an edge, it has to be differentiable.

The second term is even easier to interpret, some force wants to move us uphill or downhill in future, so simply counteract it by a new force.

If we let a counteracting forces vector act on the particles, how does this influence ? We again calculate the counteracting acceleration like we did above. Now remember that Jdescribes the slope of C with respect to any changes of the input at the present point on the multiple constraint landscapes. We get the same expression as the second term of the equation above. What we are finally looking for is an acceleration solving

We can then easily determine the counteracting force by multiplying with the mass of the particles. Are all such accelerations suitable for our problem? Let´s recall our landscape picture. When we apply a counteracting force or acceleration, this force only needs to act in the direction of the slope because we only want to counteract moving uphill and downhill. Therefore the acceleration can be represented by a coefficient times the slope of the function. There are m such constraint landscapes, so we are looking for an m-dimensional vector consisting of these coefficients that multiplied by the transpose of the slope matrix J result in the acceleration:

There is one problem left that we have to solve: Due to numerical inaccuracies the c_i will slowly go off zero despite our efforts. We have to add a spring-like mechanism to prevent these errors from accumulating. With positive c_i, we accelerate downwards and with negative c_i upwards. To prevent oscillations we add a dampening factor:

,

where k_s is the spring constant and k_dthe dampening constant.

Rearranging gives:

Now we have to solve this equation for . Normally a constraint only applies to a small number of particles (e. g. two in the case of a fixed distance), all other particle-constraint entries are zero, therefore Jis a very sparse matrix. This speeds up matrix multiplications and solving the linear equation system. After is found, it is multiplied with J^Tto get , then this force vector is added to the accumulated force vector of the particles.

Based on Siggraph 95, Course 34, Physical Based Modelling