aselab 東工大長谷川(晶)研 Physics Engine

東工大長谷川(晶)研
aselab
PHYSICS ENGINE
Rigid body simulator introduction
Continuum simulation
The Documents
http://haselab.net/class/vr
user：basic passwd：prog
Introduction for rigid body
physics engine
aselab 東工大長谷川(晶)研
What is physics engine ?
 Usual

modeling of dynamics and simulation
Manual formulation
(ml  I )  mg sin 


Find this equation with
mysterious algorithm in human brain.
1
Lagrangian also requires
(ml  I ) 2  mg cos 
2
some human intelligence
to find formulation of energy
  L  L


 (ml  I )  mg sin   0
 
t    
L  K U 
simulate this equation by ODE Solver
aselab 東工大長谷川(晶)研
What is physics engine ?

Physics engine
 Forward dynamics is automatically computed from the
configuration of the scene with formal algorithm.
fc


mv  f c  mg
  rc  f c
Iω
Other than constraint force fc ,
Newton and Euler’s equation
requires external force and it’s
application point only.
So, physics engine calculate constraint force fc with formal
algorithm.
⇒ Human intelligence is no longer required.
aselab 東工大長谷川(晶)研
Rigid body motion
v:velocity
w:angular velocity
Equation of motion
m: mass
I: inertial tensor
f: external force
mv  f
v(t  t )  v(t )  ft / m
 τ
Iω
I(t  t )ω(t  t )  I (t )ω(t )  τt
Velocity and angular velocity are constant
in case f  0, τ  0
5
aselab 東工大長谷川(晶)研
Rigid body motion simulation
 Equation
of motion
𝑓 = 𝑚𝑣
𝑁 = 𝑑(𝐼𝜔)/𝑑𝑡
 force(𝑓, 𝑁)
acts on constraints (contacts and joints)
 It is necessary to find constraint forces.
Force f
aselab 東工大長谷川(晶)研
Procedure of physics simulation
1.
Find contact points
2.
Calculate the contact forces
3.
Update the positions and velocities
Talk about this in detail
next lesson
aselab 東工大長谷川(晶)研
Taxonomy of physics simulator

Computation Methods for constraint forces are not limited to Spring and Damper.
Multi Body Dynamics simulator
Analytical method
Impulse
method
Full coordinate method
Hahn 88
Mirtich 96
Jean
Jacques
Moreau
method
LCP based formularization
Acceleration
Stewart 96
ODE

Reduced
coordinate
method
Moore &
Williams 89
Springhead1
Velocity base
base
Baraff 89
penalty method
Armstrong 79
Featherstone 83
Springhead2
OpenTissue
Velocity based LCP is widely used.

Havok Novodex ODE …
reference：
Kenny Erleben: Multibody Dynamics Animation Lecture 12
Department of Computer Science Copenhagen University
8
aselab 東工大長谷川(晶)研
Rigid body motion equation
 Two
objects＋Restraint＋External Force
Constraint
mj I j
vj wj
ri
Inertia
M
Velocity and
Angular velocity
rj
mi I i
vi wi
 Many
rigid bodies＋constraints
＋ external forces
u
force
fc
Other external
force, gravity
fe
aselab 東工大長谷川(晶)研
Constraint Examples
 Joint(Hinge) (Contacts are mentioned later)
:relative velocity/angular velocity at the joint point
w
l
j
:constraint forces and torques on the joint point
w 1 , l1
i
e1
w 4, l4
w 6 , l6
w 3, l3
e3
e2
w2,, l2
w 5 , l5
Constraint coordinate system


velocity
Only rotate the hinge axis e3 :
Rotate hinge axis e3 freely = torque not work：
force
velocity
force
10
aselab 東工大長谷川(晶)研
Embed constraints into equation of motion

Express equation of motion on the constraint coordinate.
vi ,wi: velocity or angular velocity of object i in world coordinate.
→ use w,l
w: relative velocity and angular velocity of j from i in constraint coordinate
vj ,wj
rj
w
ri
Jrow1
vi ,wi
u
Matrix expression of
cross product
Jrow4
u
express w2 w3, w5 w6, with the same way
11
aselab 東工大長谷川(晶)研
Embed constraints into equation of motion
vi ,wi: velocity, angular velocity of object i
vj , wj: velocity, angular velocity of object j
vj ,wj
rj
w
ri
w velocity/angular
u velocity/angular
vi ,wi velocity in constraint
of rigid body coordinate
coordinate
fj ,tj: constraint forces given from i to j in world coordinate
lj : constraint forces given from i to j in constraint coordinate
f j , tj
rj
l
ri
f i , ti
fc force/torque
in world coordinate
l: force/torque in constraint coordinate
12
aselab 東工大長谷川(晶)研
Embed constraints

Transform equation of motion to consider constraint easily

Ally with constraint condition (hinge’s case)
Constraint：
Substitute to equation of motion：
Solve the simultaneous equations
13
aselab 東工大長谷川(晶)研
Simulation Procedure

Calculate constraint force l

Update rigid body velocity

Update rigid body position
posture of all
rigid bodies
position
orientation
quaternion
A matrix transforming
velocity into
Quaternion differential
14
aselab 東工大長谷川(晶)研
Contact Constraint
 Contact
relative velocity at joint point
contact force on contact point


Do not penetrate each other ＝
velocity
l is not zero when w is zero or l is zero when w is not zero:
force
Static friction or kinetic friction：
velocity
force

Torque is 0：
velocity
force
15
aselab 東工大長谷川(晶)研
Embed Constraint (for contact)

Transform motion equation to consider restraint easily

Ally with constraint condition
velocity
force
velocity
force
At first, solve this case and find li.
If li does not satisfy wi = 0 then
fix li to a constant and solve wi
There are two variables in an equation at a
glance, but actually solve one variable with
consider the other variable as a constant.
16
aselab 東工大長谷川(晶)研
Summary
 Velocity

Formulate constraints into velocity based LCP and find
constraint force. Then do the stepping of the simulation




based LCP is widely used.
Transform equation of motion of solids with Jacobian J
→ Equation of motion by w and l.
Formulate constraints by equations of w and l and substitute it into
equation of motion. → The equations come to LCP
Solve LCP and find l update velocity and position.
Constraints


Joints are lines, contacts are polylines.
It looks there are two variables.
But one of the variables comes to constant.
17
aselab 東工大長谷川(晶)研
Computation
Simultaneous
Approximate
equations solution
solution for LCP
Error-correction
method
18
aselab 東工大長谷川(晶)研
How to Solve Simultaneous Equation
 Restraint
included motion equation
two variable in one equation
+complementary condition of wi and l i
(one is a variable and the other is a constant )
LCP : Linear Complementary Problem
 Solve


LCP
Pivoting method ・・・strict slow speed
Iterative method ・・・high speed, precision depends on
iteration number of times

Matrix-splitting method
Jacobi method
 Gauss–Seidel method, SOR method


Newton's method
19
aselab 東工大長谷川(晶)研
Gauss–Seidel method / SOR method
 Matrix-splitting method
A
method to solve large simultaneous equations
approximately with small computation amount.
If recurrence relation
is converged,
So
This recurrence relation converge if
A is positive symmetry matrix.
Long narrow objects or objects with small
inertia converge slow.


Choose D as a simple matrix such as diagonal matrix
Jacobi method choose D as diagonal entries of A and choose F as the rest.
Gauss-Seidel method is improved version of Jacobi method.
20
aselab 東工大長谷川(晶)研
Gauss–Seidel method / SOR method

Gauss–Seidel method
Jacobi method：update l at once
Gauss–Seidel method ：update l one row at one time
renewed
under renewing
not-renewed

SOR method
 Just
accelerate Gauss–Seidel method slightly
Gauss–Seidel method
1.6times acceleration
 Converge quickly
21
aselab 東工大長谷川(晶)研
LCP and Gauss–Seidel method
 Gauss–Seidel
method renew l one row at one time
renewed
under renewing
not-renewed
 Can
Calculate Contact force perfectly
Do not invade each other＝ l is not zero when ω is zero or l is zero
when ω is not zero:

if l1 is renewed to a negative number, then let l10

Static friction or kinetic friction：
Use renewed l1 to check “ml1 <l2< ml1”,if protrude, clip with l1m l1
22
aselab 東工大長谷川(晶)研
Simulation examples
23
Continuum simulation
aselab 東工大長谷川(晶)研
Continuum simulation

Mass point
dv(t )
f (t )  m
dt
v

t
Coil voltage
V
di (t )
V (t )  L
dt
t

Continuum
Simulation is needed not only in time but also in space dimension.
 wave
 2 h ( x, t ) 2  2 h ( x, t )
c
2
x
x
t 2
 solid
t
(i, j , k , l  x, y, z )
Discretize Hooke‘s law  ij  Eijkl kl and equation of motion
f  Ku
  Ku  Bu
Mu
 fluid
  
  
 x
y

z 
aselab 東工大長谷川(晶)研
Continuum simulation

Continuum simulation



continuous →deal with limited representative points not infinity
many methods for taking representative points
Grid points in Euler coordinate system are
representative points

Common for fluid and wave surface simulation
 Grid do not move
 Mass come from surroundings and go out

Grid points in Lagrange coordinate system
are representative points




Finite element method
Grids move with mass
Grids collapse or turn over under large deformation
Particle method： representative points are particles

Express mass with freely moving particles
 Large deformation is ok because of no assumption at beginnings
 Collision detection is important
aselab 東工大長谷川(晶)研
Finite Element Method
Soft
objects（meat, jelly, human body）
aselab 東工大長谷川(晶)研
Soft objects
cloth
Mass
point＋ Spring-damper model
Omit parts of stress
aselab 東工大長谷川(晶)研
Fluid
Liquid
Gas
aselab 東工大長谷川(晶)研
Model what?
 Wave
Particles
aselab 東工大長谷川(晶)研
Model what?
 Vortex
field of fluid
aselab 東工大長谷川(晶)研
Speed and stability
 Rigid
body stability is more important than speed
 The
number of rigid bodies is
almost enough
 Any simulation oscillates
 Guarantee control stability
Havok FX example
PD control oscillation example
Wrotek et al. 2006
32
aselab 東工大長谷川(晶)研
Speed and stability

The speed of fluid is not fast enough

2D wave front
→ sufficient real-time
3D mesh
→ more difficult
 3D particles
→ can not deal with many
 Combination is needed possibly


A Fluid Resistance Map Method for Real-time Haptic Interaction with Fluids
[Y. Dobashi et al. 2006]
3D-fluid simulation
圧力場・流速 (FRM :Fluid Resistance Map)
Prior calculation
Paddle speed
Water speed
R
Force feedback

Real-time wave simulation
Real-time
Efficient Simulation of Large Bodies of Water by Coupling Two- and Three-Dimensional Techniques
[G. Irving et al. 2006]
33
aselab 東工大長谷川(晶)研
High speed
 Soft
objects which can be simulated are numerous
 Record
and reproduce the deformation caused by impulse
force
 Modal
analysis based object deformation

Download Report