Virtual Worlds - A PARTICLE SYSTEM WITH CONSTRAINTS

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:

illustration not visible in this excerpt

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:

illustration not visible in this excerpt

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 [illustration not visible in this excerpt] 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:

illustration not visible in this excerpt

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:

illustration not visible in this excerpt

In the taylor expansion of the state vector x at the time[illustration not visible in this excerpt]

illustration not visible in this excerpt

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 [illustration not visible in this excerpt]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:

illustration not visible in this excerpt

b) Evaluate f at the midpoint of the step:

illustration not visible in this excerpt

c) Take the step using the midpoint value

illustration not visible in this excerpt

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

illustration not visible in this excerpt

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

VISCOUS DRAG

illustration not visible in this excerpt

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

SPRING FORCE

illustration not visible in this excerpt

where F a and F b are the forces on the respective particles, l = a - b, r is the rest length of the

spring[illustration not visible in this excerpt]is the derivative with respect to time, k s is 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 3 n components.

We do the same thing with the velocities, getting the derivative of our position with respect to time, that is[illustration not visible in this excerpt] We define the vector-to-vector-function C simply as the vector of all the outputs of the constraint functions:

illustration not visible in this excerpt

We assume that the initial positions are legal with all contraints, so[illustration not visible in this excerpt] The velocities are legal, too, at least for the moment, so C does not change with respect to time, that means [illustration not visible in this excerpt] In order to keep these two conditions valid from one time step to the next we have to keep [illustration not visible in this excerpt] 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 [illustration not visible in this excerpt]?

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

illustration not visible in this excerpt

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:

illustration not visible in this excerpt

where[illustration not visible in this excerpt] that is the matrix of gradients of C derived with respect to time, and[illustration not visible in this excerpt]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[illustration not visible in this excerpt] act on the particles, how does this influence[illustration not visible in this excerpt] ? We again calculate the counteracting acceleration [illustration not visible in this excerpt]like we did above. Now remember that J describes 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[illustration not visible in this excerpt] solving

illustration not visible in this excerpt

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[illustration not visible in this excerpt] consisting of these coefficients that multiplied by the transpose of the slope matrix J result in the acceleration:

illustration not visible in this excerpt

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:

illustration not visible in this excerpt

where k s is the spring constant and k d the dampening constant. Rearranging gives:

illustration not visible in this excerpt

Now we have to solve this equation for [illustration not visible in this excerpt] 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 J is a very sparse matrix. This speeds up matrix multiplications and solving the linear equation system. After[illustration not visible in this excerpt] is found, it is multiplied with[illustration not visible in this excerpt]to get[illustration not visible in this excerpt] then this force vector is added to the accumulated force vector of the particles.

Based on Siggraph 95, Course 34, Physical Based Modelling