Particle dynamics Particle overview Particle system Forces - - PowerPoint PPT Presentation

Particle dynamics

• Particle overview
• Particle system
• Forces
• Constraints
• Second order motion analysis

Particle system

• Particles are objects that have mass, position,

and velocity, but without spatial extent

• Particles are the easiest objects to simulate but

they can be made to exhibit a wide range of

• bjects

A Newtonian particle

• First order motion is sufficient, if
• a particle state only contains position
• no inertia
• particles are extremely light
• Most likely particles have inertia and are

affected by gravity and other forces

• This puts us in the realm of second order

motion

Second-order ODE

f = ma

What is the differential equation that describes the behavior of a mass point? What does f depend on?

¨ x(t) = f(x(t), ˙ x(t)) m

Second-order ODE

Add a new variable, v(t), to get a pair of coupled first order equations

This is not a first oder ODE because it has second derivatives

{ ˙

x = v ˙ v = f/m

¨ x(t) = f(x(t), ˙ x(t)) m = f(x, ˙ x)

Phase space

• x

v

• =

⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ x1 x2 x3 v1 v2 v3 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦

Concatenate position and velocity to form a 6-vector: position in phase space First order differential equation: velocity in the phase space

˙ x ˙ v

• =

v

f m

Quiz

• A mass point attached to a spring obeys Hooke’s

Law: f = -k (x - x0). What is the ODE that describes this motion?

• Particle overview
• Particle system
• Forces
• Constraints
• Second order motion analysis

Particle structure

x v f m

Particle

position velocity force accumulator mass

a point in the phase space

Solver interface

system solver solver interface

x v f m

particle GetDim

6

Get/Set State

x v

Deriv Eval

v

f m

system

Particle system structure

n time ...

x1 v1 f1 m1 x2 v2 f2 m2 xn vn fn mn

particles

system

Particle system structure

solver solver interface

GetDim

6n

Get/Set State

x1 v1 x2 v2 xn vn

. . .

Deriv Eval

vn

fn mn

v1

f1 m1

v2

f2 m2

. . .

particles time n

Deriv Eval

Clear forces: loop over particles, zero force accumulator Calculate forces: sum all forces into accumulator Gather: loop over particles, copy v and f/m into destination array

• Particle overview
• Particle system
• Forces
• Constraints
• Second order motion analysis

Forces

• Constant
• gravity
• Position dependent
• force fields, springs
• Velocity dependent
• drag

• system

particles n time

Particle systems with forces

...

x1 v1 f1 m1 x2 v2 f2 m2 xn vn fn mn

F1 F2 Fm ... forces

Force structure

• Unlike particles, forces are heterogeneous

(type-dependent)

• Each force object “knows”
• which particles it influences
• how much contribution it adds to the force

accumulator

system particles n time

Particle systems with forces

...

x1 v1 f1 m1 x2 v2 f2 m2 xn vn fn mn

forces F1 F2 Fm ...

F

Gravity

p sys apply_fun G particle system

. . .

x1 v1 f1 m1 x2 v2 f2 m2 xn vn fn mn

p->f += p->m*F->G

Unary force: f = mG Exerting a constant force on each particle

F

Viscous drag

k sys particle system p

. . .

x1 v1 f1 m1 x2 v2 f2 m2 xn vn fn mn

apply_fun p->f += p->v*F->k

At very low speeds for small particles, air resistance is approximately:

fdrag = −kdragv

Attraction

Act on any or all pairs of particles, depending on their positions xp xq l

fq = −fp fp = −k mpmq |l|2 l |l| l = xp − xq

Attraction

p sys apply_fun k particle system

F

xp vp fp mp xq vq fq mq

fp = −k mpmq |l|2 l |l|

Damped spring

fp = −

• ks(|l| − r) + kd

˙ l · l |l|

• l

|l|

r | l | xp xq

l = xp − xq fq = −fp

Damped spring

p sys apply_fun ks particle system

F

xp vp fp mp xq vq fq mq

kd r

fp = −

• ks(|l| − r) + kd

˙ l · l |l|

• l

|l|

Quiz

For an ideal spring, what is the force it applies to two particles, p and q, attached to it. Write down the pseudo code for its “apply_fun”.

p sys apply_fun ks particle system

F

r

xp vp fp mp xq vq fq mq

Deriv Eval

...

x1 v1 f1 m1 x2 v2 f2 m2 xn vn fn mn

• 1. Clear force accumulators

F1 F2 Fm ...

• 2. Invoke apply_force

functions

• 3. Return derivatives to solver

˙ x ˙ v

• =

v

f m

ODE solver

Euler’s method:

x(t0 + h) = x(t0) + hf(x, t) xt+1 = xt + h ˙ xt vt+1 = vt + h ˙ vt

2. 4.

Get/Set State

Euler step

system

GetDim Deriv Eval

solver solver interface

xt+1 = xt + h ˙ xt vt+1 = vt + h ˙ vt

3. time

• 5. Advance

time Deriv Eval

1.

. . .

vn

fn mn

v1

f1 m1

v2

f2 m2

x1 v1 x2 v2 xn vn

. . .

particles

• Particle overview
• Particle system
• Forces
• Constraints
• Second order motion analysis

Particle Interaction

• We will revisit collision

when we talk about rigid body simulation

• For now, just simple

point-plane collisions

Collision detection

vT vN v x N

Normal and tangential components vN = (N · v)N vT = v − vN

Collision detection

vT vN v

p x

N

Particle is on the legal side if Particle is heading in if

Particle is within of the wall if

• (x − p) · N ≥ 0

(x − p) · N < v · N < 0

Collision response

v vT vN

Before collision

v′ vT −krvN

After collision coefficient of restitution:

0 ≤ kr < 1

v′ = vT − krvN

Contact

p

N

Conditions for resting contact:

f

fN fT

If a particle is pushed into the contact plane a contact force fc is exerted to cancel the normal component of f

x

• 1. particle is on the collision surface

v

• 2. zero normal velocity

• Particle overview
• Particle system
• Forces
• Constraints
• Second order motion analysis

Linear analysis

• Linearly approximate acceleration
• Split up analysis into different cases
• constant acceleration
• linear acceleration

d dt  x v

• =

 I −K −D  x v

• + a0

d dt  x v

• = f(

 x v

• , t) = A

 x v

• + a0

Constant acceleration

• Solution is
• v(t) only needs 1st order accuracy, but x(t)

demands 2nd order accuracy

v(t) = v0 + a0t x(t) = x0 + v0t + 1 2a0t2

Linear acceleration

• When K (or D) dominates ODE, what type of

motion does it correspond to?

• Need to compute the eigenvalues of A

d dt

• x

v

• =
• I

−K −D x v

• = A
• x

v

Linear acceleration

Assume is an eigenvalue of A, is the corresponding eigenvector

α

• u1

u2

• I

−K −D u1 u2

• = α
• u1

u2

• The eigenvector of A has the form
• u

αu

• Assuming D is linear combination of K and I (Rayleigh

damping) That means K and D have the same eigenvectors

Linear acceleration

Now assume u is an eigenvector for both K and D For any u, if is an eigenvector of A, the following must be true

• u

αu

• I

−K −D u αu

• = α
• u

αu

• −λku − αλdu = α2u

α = −1 2λd ±

• (1

2λd)2 − λk

Eigenvalue approximation

• If D dominates
• exponential decay
• If K dominates
• oscillation

α ≈ −λd, 0 α ≈ ± √ −1

• λk

Analysis

• Constant acceleration (e.g. gravity)
• demands 2nd order accuracy for position
• Position dependence (e.g. spring force)
• demands stability, oscillatory motion
• looks at imaginary axis
• Velocity dependence (e.g. damping)
• demands stability, exponential decay
• looks at negative real axis

Explicit methods

• First-order explicit Euler method
• constant acceleration:
• position dependence:
• velocity dependence:
• RK3 and RK4
• constant acceleration:
• position dependence:
• velocity dependence:

bad (1st order) very bad (unstable)

• k (conditionally stable)

great (high order)

• k (conditionally stable)
• k (conditionally stable)

Implicit methods

• Implicit Euler method
• constant acceleration:
• position dependence:
• velocity dependence:
• Trapezoidal rule
• constant acceleration:
• position dependence:
• velocity dependence:

bad (1st order)

• k (stable but damped)

great (monotone) great (2nd order) great (stable and no damp) good (stable, not monotone)