Particle system dynamics
Jiří Chmelík, Marek Trtík
PA199
► Motivation
► Motion of a single particle: Equations of motion
► Use of an ODE solver
► Motion of many particles
► Forces
► Gravity, drag, spring, local interaction
► Collision: particle vs. plane
► Detection, response, simple friction
Outline
Motivation
3https://experiments.withgoogle.com/fluid-particles
https://en.wikipedia.org
/wiki/Particle_system
https://github.com/Lakshitha
Madushan/Unity-Particle-
System
► Particle = an abstract object with these properties:
► No spatial extent - it is just a point in 3D space
► Velocity
► Respond to forces (e.g., gravity)
► Mass - resistance to changes in motion state
► Particle in math: 𝒫 = (𝒙, 𝒗, 𝑭, 𝑚).
► Particle in C++:
struct Particle {
Vector3 position;
Vector3 velocity;
Vector3 force;
float mass;
};
Particle definition
4
► Motion of a particle 𝒫 in space is given by a function of time:
► 𝒫(𝑡)= 𝒙, 𝒗, 𝑭, 𝑚 𝑡 = (𝒙 𝑡 , 𝒗 𝑡 , 𝑭, 𝑚)
► 𝑚 is constant (not dependent on time).
► 𝑭 is total external force (not updated by the particle system).
► To compute 𝒫(𝑡) we need to know how it changes in time.
=> We need to compute ሶ𝒫 𝑡 = ( ሶ𝒙 𝑡 , ሶ𝒗 𝑡 ).
► Newton’s second law of motion: 𝑭 = 𝑚𝒂
► Important relations: 𝒗 = ሶ𝒙 =
𝑑𝒙
𝑑𝑡
, 𝒂 = ሶ𝒗 =
𝑑𝒗
𝑑𝑡
.
►So, 𝒫(𝑡) is a solution of Newton’s equations of motion:
ሶ𝒙 𝑡 = 𝒗 𝑡 , ሶ𝒗(𝑡) = 𝒂 =
𝑭
𝑚
.
Particle equations of motion
5
► There is 6 ordinary differential equations (ODE) of the 1st order in the
Newton’s equations of a single particle.
► 𝒙 𝑡 and 𝒗 𝑡 are 3D vector functions.
► In general, a system of 𝑛 1st order ODEs has the form:
ሶ𝒚 = 𝑭(𝒚, 𝑡)
where 𝒚(𝒕) = 𝑦0 𝑡 , … , 𝑦 𝑛−1 𝑡
⊤
and
𝑭(𝒚, 𝒕) = 𝐹0 𝑦0 𝑡 , … , 𝑦 𝑛−1 𝑡 , t , … , 𝐹𝑛−1 𝑦0 𝑡 , … , 𝑦 𝑛−1 𝑡 , t
⊤
.
Therefore, we have a system:
ሶ𝑦0 = 𝐹0 𝑦0, … , 𝑦 𝑛−1, 𝑡 , … , ሶ𝑦 𝑛−1 = 𝐹𝑛−1 𝑦0, … , 𝑦 𝑛−1, 𝑡
►At each simulation time 𝑡0 we know 𝒙 𝑡0 = 𝑿0 and 𝒗 𝑡0 = 𝑽0.
►Therefore, we solve the initial value problem of 1st order ODEs:
ሶ𝒚 = 𝑭(𝒚, 𝑡), 𝒚 𝑡0 = 𝒚0
Solving equations of motion
6
► We are given a black-box function ODE solving the initial value
problem of 1st order ODEs ሶ𝒚 = 𝑭(𝒚, 𝑡), 𝒚 𝑡0 = 𝒚0 :
using F_y_t = std::function<float(std::vector<float> const&,float)>;
void ODE(
std::vector<float> const& y0, // 𝑿0, 𝑽0 of particle(s)
std::vector<F_y_t> const& Fyt, // ሶ𝒙, ሶ𝒗 of particle(s), i.e. 𝒗, 𝑭/𝑚
float& t, // current time (to be updated)
float const dt, // time step
std::vector<float>& y // integrated 𝒙, 𝒗 of particle(s)
);
NOTE: Implementation of ODE is the topic of next lecture.
Solving equations of motion
7
void getState(Particle const& p, std::vector<float>& y0) {
y0.push_back(p.position.x);
y0.push_back(p.position.y);
y0.push_back(p.position.z);
y0.push_back(p.velocity.x);
y0.push_back(p.velocity.y);
y0.push_back(p.velocity.z);
}
Building initial state for ODE
8
void getDerivative(Particle const& p, std::vector<F_y_t>& Fyt) {
Fyt.push_back([&p](std::vector<float> const&,float){ return p.velocity.x; });
Fyt.push_back([&p](std::vector<float> const&,float){ return p.velocity.y; });
Fyt.push_back([&p](std::vector<float> const&,float){ return p.velocity.z; });
Fyt.push_back([&p](std::vector<float> const&,float){ return p.force.x/p.mass; });
Fyt.push_back([&p](std::vector<float> const&,float){ return p.force.y/p.mass; });
Fyt.push_back([&p](std::vector<float> const&,float){ return p.force.z/p.mass; });
}
►Observation: Parameters of lambda functions are not used.
►Our functions 𝑭 𝒚, 𝑡 are simple; ODE solver handles general case.
Building derivatives for ODE
9
void doSimulationStep(Particle& p, float& t, float const dt) {
UpdateForce(p,t,dt); // Applies external forces and impulses.
std::vector<float> y0, y;
std::vector<F_y_t> Fyt;
getState(p, y0);
getDerivative(p, Fyt);
ODE(y0, Fyt, t, dt, y); // Computes y and updates t (t += dt).
setState(p, y.begin());
}
Simulation step for single particle
10
void setState(Particle& p, std::vector<float>::const_iterator& it) {
p.position.x = *it; ++it;
p.position.y = *it; ++it;
p.position.z = *it; ++it;
p.velocity.x = *it; ++it;
p.velocity.y = *it; ++it;
p.velocity.z = *it; ++it;
}
Saving ODE results
11
Data flow in simulation step
12
𝑥 𝑦 𝑧 𝑥 𝑦 𝑧
𝑓𝑥 𝑓𝑦 𝑓𝑧 𝑓𝑥 𝑓𝑦 𝑓𝑧
𝑥 𝑦 𝑧 𝑥 𝑦 𝑧
y0
Fyt
y
ODE
p.position p.velocity
p.position p.velocity
p.force / p.massp.velocity
► It is a system consisting of 𝑛 patricles.
► Particle system in math:
𝒫 𝑛 = 𝒫0, 𝒫1, … , 𝒫𝑛−1 =
[ 𝒙0, 𝒗0, 𝑭0, 𝑚0 , 𝒙1, 𝒗1, 𝑭1, 𝑚1 , … , 𝒙 𝑛−1, 𝒗 𝑛−1, 𝑭 𝑛−1, 𝑚 𝑛−1 ].
► Particle system in C++:
using ParticleSystem = std::vector<Particle>;
Particle system
13
void getState(ParticleSystem const& ps, std::vector<float>& y0) {
for (Particle const& p : ps) getState(p,y0);
}
void getDerivative(ParticleSystem const& ps, std::vector<F_y_t>& Fyt) {
for (Particle const& p : ps) getDerivatives(p, Fyt);
}
void setState(ParticleSystem& ps, std::vector<float>::const_iterator& it) {
for (Particle& p : ps) setState(p, it);
}
ODE helper functions
14
void doSimulationStep(ParticleSystem& ps, float& t, float const dt) {
UpdateForce(ps,t,dt); // Applies external forces and impulses.
std::vector<float> y0, y;
std::vector<F_y_t> Fyt;
getState(ps, y0);
getDerivative(ps, Fyt);
ODE(y0, Fyt, t, dt, y); // Computes y and updates t (t += dt).
setState(ps, y.begin());
}
Simulation step for whole system
15
Data flow in simulation step
16
NOTE: For 𝒫 𝑛 we have a system of 6𝑛 equations.
𝑥 𝑦 𝑧 𝑥 𝑦 𝑧
𝑓𝑥 𝑓𝑦 𝑓𝑧 𝑓𝑥 𝑓𝑦 𝑓𝑧
𝑥 𝑦 𝑧 𝑥 𝑦 𝑧
y0
Fyt
y
ODE
𝑥 𝑦 𝑧 𝑥 𝑦 𝑧
𝑓𝑥 𝑓𝑦 𝑓𝑧 𝑓𝑥 𝑓𝑦 𝑓𝑧
𝑥 𝑦 𝑧 𝑥 𝑦 𝑧
𝑥 𝑦 𝑧 𝑥 𝑦 𝑧
𝑓𝑥 𝑓𝑦 𝑓𝑧 𝑓𝑥 𝑓𝑦 𝑓𝑧
𝑥 𝑦 𝑧 𝑥 𝑦 𝑧
ps[0] ps[1] ps[2]
…
void UpdateForce(ParticleSystem& ps, float const t, float const dt) {
clearForce(ps);
applyForce(ps,t,dt); // Add all forces and impulses to all particles.
}
void clearForce(ParticleSystem& ps) {
for (Particle& p : ps) p.force = Vector3(0,0,0);
}
►Next we discuss what forces we can add to particles inside the
function applyForce().
Forces
17
► Homogenous field:
► For each particle we add the force vector 𝑭 = 𝑚𝒈 where
► 𝑚 is the mass of the particle.
► 𝒈 is a constant vector, e.g., 𝒈 = Vector3(0,0, −10).
► Radial field:
► There is a center of gravity 𝑺 of mass 𝑀 (it can be one of the particles).
► For each particle we add the force vector
𝑭 = 𝐺
𝑀𝑚
𝑺−𝒙 2
𝑺−𝒙
𝑺−𝒙
= 𝐺
𝑀𝑚
𝑺−𝒙 3 (𝑺 − 𝒙), where
► 𝐺 is the gravitational constant.
► 𝑚 is the mass of the particle.
► 𝒙 is the position of the particle.
► We can handle cases when 𝑺 − 𝒙 is small by not applying the force.
Gravity
18
► A force of the environment making a particle decrease its velocity
relative to the environment.
► A drag force can also enhance numerical stability of simulation.
► For each particle we add the force vector 𝑭 = 𝑘 𝑑(𝑽 − 𝒗), where
► 𝑘 𝑑 is the coefficient of drag.
► 𝑽 is the velocity of the environment (often 𝑽 = 𝟎).
► 𝒗 is the velocity of the particle.
Viscous Drag
19
► It is a force between two particles 𝒫𝑖 and 𝒫𝑗 given by Hook’s law:
𝑭𝑖 = − 𝑘 𝑠 𝒅 − 𝑑0 + 𝑘 𝑑
ሶ𝒅 ⋅
𝒅
𝒅
𝒅
𝒅
𝑭𝑗 = −𝑭𝑖 (3rd Newton’s law – action and reaction)
where
► 𝑘 𝑠 is the spring constant.
► 𝑘 𝑑 is the damping constant.
► 𝒅 = 𝒙𝑖 − 𝒙𝑗 is the distance vector between the particles.
► 𝑑0 is the rest length between the particles.
► ሶ𝒅 = 𝒗𝑖 − 𝒗𝑗 is the relative velocity between the particles.
Spring
20
► Particles start to interact when they come close.
► Particles stop to interact when they move apart.
► Example: Particle-based fluid simulation.
► Computationally expensive task:
► 𝒪(𝑛2
) – all pairs of particles are checked.
► Space partitioning methods (e.g., octree)
are essential for performance.
Local interaction
21
https://experiments.withgoogle.com/fluid-particles
► We often want particles to collide with the ground or a wall. These
boundaries can be approximated by planes.
► The process consists of two parts:
► Detection of a collision.
► Response to the collision.
Collision: particle vs. plane
22
https://github.com/LakshithaMadushan/Unity-
Particle-System
► Let us consider a particle 𝒫 = (𝒙, 𝒗, 𝑭, 𝑚).
► The plane is represented by the equation 𝑵 ⋅ 𝑿 − 𝑷 = 0, where
► 𝑵 is the unit normal vector pointing “outside” (above the ground).
► 𝑷 is some point in the plane.
► 𝑿 is a tested point.
► The particle collides with the plane only if 𝑵 ⋅ 𝒙 − 𝑷 ≤ 0.
► Only in that case we proceed to the collision response.
Collision detection
23
► If the particle increases the penetration with the plane, i.e., when 𝑵 ⋅
𝒗 < 0, then we change the component of 𝒗 orthogonal to the plane:
► The component of 𝒗 orthogonal to the plane is 𝒗⊥
= 𝑵 ⋅ 𝒗 𝑵.
► The velocity change is then is Δ𝒗 = − 1 + 𝑟 𝒗⊥
= − 1 + 𝑟 𝑵 ⋅ 𝒗 𝑵, where
► 𝑟 ∈ 0,1 is the coefficient of restitution.
► We update 𝒗 to be 𝒗 + Δ𝒗.
► NOTE: Formally, we apply an impulse 𝑰 = 𝑚Δ𝒗 to the particle.
► If 𝑵 ⋅ 𝑭 < 0, then we cancel the component of 𝑭 orthogonal to the
plane:
► We compute Δ𝑭 = −𝑭⊥, where 𝑭⊥ = 𝑵 ⋅ 𝑭 𝑵.
► We update 𝑭 to be 𝑭 + Δ𝑭.
► NOTE: This step should be applied after all external forces (gravity, etc.)
were added to the 𝑭 field of the particle.
Collision response
24
► We build a simplified friction model for particle system:
► We do not distinguish static and dynamic friction.
► We ignore variable changes caused by interactions with other particles.
► If 𝑵 ⋅ 𝑭 < 0, then a friction force 𝑭 𝑓 is acting on the particle:
► |𝑭 𝑓| is proportional to |𝑭⊥
|, where 𝑭⊥
= 𝑵 ⋅ 𝑭 𝑵.
► The direction of 𝑭 𝑓 is opposite to the component 𝒗∥
of 𝒗 parallel with the
plane, where 𝒗∥
=
𝑵×𝒗×𝑵
𝑵×𝒗×𝑵
.
► Therefore, we define the friction force as 𝑭 𝑓 = 𝑘 𝑓 𝑵 ⋅ 𝑭 𝒗∥, where
► 𝑘 𝑓 is a friction coefficient.
► Note: We should apply the friction before the collision response.
Simple friction
25
► We defined particle and particle system.
► We learned Newton’s equations of motion for a particle, i.e., a system
of 1st order ODEs.
► We learned how to use ODE solver for the simulation.
► We learned several kinds of forces which we can apply to particles.
► We know how to compute and respond to collision of a particle with a
plane, including application of a friction force.
Summary
26
► [1] Andrew Witkin; Physically Based Modeling: Principles and Practice
Particle System Dynamics; Robotics Institute, Carnegie Mellon
University, 1997.
References
27