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/s2 ).
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.
Frequently asked questions
What is a particle system and how does it work?
A particle system is a collection of particles, each with properties like mass, position, and velocity. These particles move based on forces acting on them, such as gravity, drag, or spring forces. They can model various physical phenomena, including dust, snow, chains, and rigid bodies, by restricting particle movement appropriately.
How is the position of a particle updated in a particle system?
The position of a particle is updated using the following formula: x(t+dt) = x(t) + v(t)dt + 0.5 * F(t) / m * dt2, where x is the position vector, v is the velocity vector, F is the force acting on the particle, m is the mass, and dt is the time step.
What is the ordinary differential equation used to describe the behavior of a particle system?
The behavior of the particle system is described by the equation: dx/dt = f(x), where x is a 6n-dimensional vector containing the positions and velocities of all particles, and f determines the derivative of each component of x.
What are some methods for numerically solving differential equations in particle systems?
Several methods exist, including Euler's method, the midpoint method, and Runge-Kutta methods. Euler's method is a simple but less accurate method. The midpoint method is more accurate. Runge-Kutta methods offer even higher accuracy but require more computation.
What are some common forces that act on particles in a particle system?
Common forces include gravity (F = m * g), viscous drag (F = -k d * v), and spring force (F a = -k s * (|l| - r) * l/|l| - k d * (l´) ).
Why are constraints important in particle systems?
Constraints allow for more complex and realistic simulations by defining conditions for legal particle positions, such as fixing particles in space or maintaining a fixed distance between them.
What is the problem with using firm springs to implement constraints?
Firm springs oscillate very fast, requiring very small time steps to be simulated correctly, which consumes a lot of computational power. Instabilities can also occur, disrupting the simulation.
How are constraints generally defined in this system?
Constraints are defined using m functions c i, one for each constraint, that must be zero when the respective constraint is met. These functions must be differentiable and take the positions of all particles as input.
What is the Jacobian matrix J and how is it used?
The Jacobian matrix J of C contains the gradients of the constraint functions c i. It describes the change of C when parameters of C are changed, and is used to determine the forces needed to counteract constraint violations.
How are numerical inaccuracies addressed in constraint implementation?
A spring-like mechanism is added to prevent errors from accumulating due to numerical inaccuracies. This mechanism applies forces to correct deviations from the constraint conditions, with a dampening factor to prevent oscillations: δ a = -J T[k s C + k d (Jq´)].
- Quote paper
- Arno Schödl (Author), 1997, Virtual Worlds - A PARTICLE SYSTEM WITH CONSTRAINTS, Munich, GRIN Verlag, https://www.hausarbeiten.de/document/96310
