Modelagem Baseada na Física
Simulação de Corpos Rígidos
César Candido Xavier
Mestrado Computação Gráfica
UFRJ - COPPE
1
Objetivo
Apresentar os principais conceitos das notas
de aula elaboradas por David Baraff no
Modelamento Baseado na Física para
simulação de movimentos de Corpos
Rígidos.
2
Roteiro
• Modelando o Movimento de uma Partícula
• Modelando o Movimento um Corpo Rígido
–
–
–
–
–
–
•
•
•
•
Velocidade Angular
Massa de um Corpo
Centro de Massa
Força e Torque
Momento Linear
Momento Angular
• Momento Inércia
Equação Movimento de um Corpo Rígido
Codificação em C++
Conclusão
Bibliografia
3
Modelando o Movimento de uma Partícula
Movimento de uma Partícula
4
Estado da Partícula
 x(t ) 

Y  
 v(t ) 
5
Dinâmica da Partícula
6
Variáveis de Estado
d
d  x(t )   v(t ) 
Y  
  

dt
dt  v(t )   F (t ) / m 
7
Múltiplas Partículas
v(t )
F (t )
F (t )
v(t )
v(t )
F (t )
F (t )
v(t )
8
Variáveis de Estado
 x1 (t )   v1 (t ) 

 

 v1 (t )   F1 (t ) / m 
d
d   


Y

dt
dt 
 

x
(
t
)
v
(
t
)
n
n

 

 v (t )   F (t ) / m 
 n   n

9
Solução EDO
10
Solução EDO
Y(t0)
len
t0
Solucionador
Y(t1)
EDO
t1
dydt
void dydt(double t, double y[ ],double ydot[ ])
11
dydt
x1 (t )
v1 (t )

Y (t ) 
xn (t )
vn (t )
d
Y (t ) 
dt
v1 (t )
F1 (t ) / m

vn (t )
Fn (t ) / m
12
Modelando um Corpo Rígido
• Variáveis de Estado do Corpo Rígido
13
Equação Movimento de um Corpo Rígido
 x1 (t )   v1 (t ) 

 

d
d ?   ? 
Y 

dt
dt Mv(t )   F (t ) 

 

 ?   ? 
14
Malha de Forças
f2
f1
f1
F (t )   f 1
f3
f1
15
Orientação
• Iremos representar a orientação de um
corpo rígido pela matriz de rotação R(t).
Os pontos são transformados de
coordenadas do corpo para coordenadas
do espaço como:
P(t)=R(t)p0 + x(t)
16
Mudança de Coordenadas
Coordenadas Espaço Corpo
Coordenadas Espaço
17
Velocidade Angular
• Representamos a velocidade angular como um
vetor, w(t), o qual codifica tanto a velocidade
quanto o eixo de giro.
– Como w(t) e R(t)
se relacionam ?
18
Interpretação Física da Matriz R(t)
 rxx

• Sabemos que: R(t )   rxy
r
 xz
ryx
ryy
ryz
rzx 

rzy 
rzz 
1
 
• Quem seria: R(t ) 0  ?
 0
 
 1   rxx 
   
R (t ) 0    rxy 
0  r 
   xz 
19
Relação entre w(t) e R(t)
• r(t) é fixo ao corpo, logo independe da translação
do corpo;

• r (t )  w(t ) b

• r (t )  w(t )  b
20
Relação entre w(t) e R(t)
• Sabemos que no tempo t a direção do eixo x é
dado pela primeira coluna da matriz de rotação.
• A derivada desta, é a taxa de mudança deste
vetor que é:
r
 xx 
 
w(t )   rxy 
r 
 xz 
• Sabemos entretanto que a  b é:
 a y bz  by a z 


  a x bz  bx a z 
 a b b a 
 x y x y 
21
Relação entre w(t) e R(t)
• Dado um vetor a podemos definir a* como sendo a
matriz:
 0


a   az
 a
 y
 az
0
ax
• Então a  b pode ser escrito:
 0


a b   az
 a
 y
 az
0
ax
ay 

 ax 
0 
a y  bx   a y bz  by az 
  

 ax  by     ax bz  bx az   a  b
0  bz   ax by  bx a y 
22
Relação entre w(t) e R(t)
• Utilizando a formulação anterior temos que:

 rxx 


 
R(t )   w(t )   rxy 

r 
 xz 

 ryx 
 
w(t )   ryy 
r 
 yz 
• Pode ser escrito como:
 rxx


R(t )  w(t )  rxy
r
 xz

ryx
ryy
ryz
 rzx  
 
w(t )   rzy  
 r 
 zz  
rzx 

rzy 
rzz 
• Como a matriz à direita é R(t) podemos escrever:

R(t )  w(t ) R(t )
23
Massa de um Corpo
• Consideremos um corpo rígido feito de um
número grande de pequenas partículas, cada
uma localizada:
ri (t )  R(t )r0 i  x(t )
• A massa total M do corpo será
N
M   mi
i 1
24
Velocidade de uma Partícula
• A velocidade é dada pela derivada da posição
da partícula, ou seja:

r i (t )  w(t ) R(t )r0i  v(t )
• Podemos reescrever como:

r i (t )  w(t ) ( R(t )r0i  x(t )  x(t ))  v(t )

r i (t )  w(t ) (ri (t )  x(t ))  v(t )

r i (t )  w(t )  (ri (t )  x(t ))  v(t )
25
Velocidade de uma Partícula
26
Centro de Massa
N
• É definido como  m r (t )
i i
i 1
M
• Quando utilizamos o centro de massa como
sistema de coordenadas do corpo, queremos
dizer que:
N
 m r (t )
i 1
i i0
M
 0
 
 0   0
 0
 
27
Centro de Massa
• Qual a localização do centro de massa em t ?
N
N
 m r (t )  m ( R(t )r (t )  x(t ))
i 1
i i
M

i 1
i
0i
M

N
N
i 1
i 1
R(t ) mi r0 i   mi x(t )
M
N
 x(t )
m
i 1
M
i
 x(t )
• Outra relação importante:
 m (r (t )  x(t ))  0
i
i
28
Equação do Movimento Corpo Rígido
 x1 (t )   v1 (t ) 

 


d
d  R(t )   w(t ) R(t ) 
Y 

dt
dt Mv(t )   F (t ) 

 

?
 ?  

29
Força e Torque
• Torque difere da força uma vez que o
torque em uma partícula depende da
localização relativa da mesma em relação
ao centro de massa x(t).
• Intuitivamente a direção do torque é o do
eixo em torno do qual o corpo gira.
• F(t) não traz informações sobre onde as
várias forças agem em um corpo, ao
contrário de τ(t), que nos dá a distribuição
de forças sobre o corpo.
30
Força e Torque
31
Força e Torque
32
Momento Linear
• O Momento Linear de uma partícula é definido
como sendo: p  mv

• O momento total P(t) é: P(t )   mi r i (t )
• Podemos escrevê-lo como:
– P(t )   (mi v(t )  mi w(t )  (ri (t )  x(t )))
– P(t )   (mi v(t )   mi w(t )  (ri (t )  x(t )))
– P(t )   mi v(t )   mi  v(t )  Mv(t )
33
Momento Linear
• O momento linear de um corpo rígido é o mesmo
se ele fosse uma partícula de massa M e
velocidade v(t).
• Obtemos diretamente:

–

v(t ) 
P(t )
M

–
P(t )  F (t )
–

F (t )
v(t ) 
M
34
Momento Angular
• É o menos intuitivo dos conceitos vistos até
aqui...
• Proporcionará equações mais simples...
• O momento angular é definido como:
l  r  p, cujo módulo é l  rpsen( ).
35
Momento Angular


r
uv
p

36
Momento Angular
• Sabemos que:

p(t )  F (.t )
• Multiplicando vetorialmente por r :
–


r  F (t )  r  p(t ) ou seja,   r  p(t )
dl d (r  p)

• Sabemos que:
dt
dt
, o qual podemos escrever como:
dl dr
dp


p

r

–
dt dt
dt
dl
dp
 (v  mv )  r 
–
dt
dt

dl
dp
 r

dt
dt
37
Momento Angular
• O momento angular total de um corpo rígido
será:
N
L  l1  l2    l N   li
i 1
dL(t )
• E o torque total será:  
d (t )
dP(t )
• Que é o equivalente de F 
d (t )
38
Tensor (Momento) de Inércia
• Em um corpo rígido as partículas mantêm as
mesmas posição relativas.
• Consideremos um corpo rígido girando com
velocidade angular em torno de um eixo. Seja K a
energia cinética deste corpo.
K
1
m1v12  m2v2 2    mN vN 2 ,
2
como vi  wri , temos :
1
1 N
2
2
2
2
K  m1 ( wr1 )  m2 ( wr2 )    mN ( wrN )   ( mi ri ) w2
2
2 i 1
N
I   mi ri
2
i 1
39
Relação entre o Tensor de Inércia e a
Velocidade Angular
• Podemos demonstrar que:
L(t )  I (t )w(t )
• Que é a forma similar de P(t )  Mv(t )
• O momento de inércia I(t) é o fator de escala entre
o momento angular L(t) e a velocidade angular
w(t).
40
Momento de Inércia
• Seja ri’ o deslocamento da i-ésima partícula em
relação a x(t) no tempo t, ou seja:
ri  ri (t )  x(t )
'
• O momento de inércia I(t) é expresso em termos
de ri’ como a seguinte matriz simétrica:
 mi rix´ riy´
 mi rix´ riz´ 
 mi (riy´2  riz´2 )


´ ´
´2
´2
´ ´
I (t )     mi riy rix
mi (rix  riz )
 mi riy riz 
´ ´
´2
´2 
  m r´r´

m
r
r
m
(
r

r
)
i iz ix
i iz iy
i
iy
ix

41
Momento de Inércia
• Usando o fato de que r r  r  r  r
I(t) como a diferença :
'T
i
'
i
´2
ix
´2
iy
´2
iz
´2
 1 0 0   mi rix

  ´ ´
´T ´
I (t )   mi ri ri  0 1 0   mi riy rix
 0 0 1   m r´r´

  i iz ix
podemos reescrever
mi rix´ riy´
mi riy´2
mi riz´ riy´
mi rix´ riz´ 

´ ´
mi riy riz 
mi riz´2 
• Tomando o produto:
 rix 
 ´ ´
' 'T
ri ri   riy rix
 r´ 
 iz 
´
riy´
 rix´2
 ´ ´
´
riz    riy rix
 r´r´
 iz ix
rix´ riy´
riy´2
riz´ riy´
rix´ riz´ 

´ ´
riy riz 
riz´2 
42
Momento de Inércia
• Seja 1 a matriz unitária 33, podemos expressar
I(t) como:
I (t )   mi ((ri ri )1 - ri ri )
'T
'
'
'T
r0 i
• Sabemos que ri (t )  R(t )r0 i  x(t ) onde
é
´
constante. Daí ri  R(t )r0 i e podemos escrever I(t)
como:
I (t )   mi ((R(t )r0 i ) ( R(t )r0 i )1 - ( R(t )r0 i )( R(t )r0 i ) )
T
T
I (t )   mi (r0 i R(t ) R(t )r0 i )1 - ( R(t )r0 i r0 i R(t ) )
T
T
T
T
43
Momento de Inércia
I (t )   mi (r0 i R(t ) R(t )r0 i )1 - ( R(t )r0 i r0 i R(t ) )
T
T
T
T
I (t )   mi (r0 i r0 i )1 - ( R(t )r0 i r0 i R(t ) )
T
T
T
I (t )   mi (( R(t )(r0Ti r0 i ) R(t )T 1 - ( R(t )r0 i r0Ti R(t )T )
I (t )  R(t )( mi ((r0 i r0 i )1 - r0 i r0 i )) R(t )
T
T
T
I body   mi ((r0i r0i )1 - r0i r0i )
T
I (t )  R(t )( I body ) R(t )
T
T
44
Inverso do Momento de Inércia
I (t )  ( R(t )( I body ) R(t ) )
1
T
I (t )  ( R(t ) ) I body R(t )
1
T
1
1
I (t )  R(t ) I body R(t )
1
1
1
1
T
45
Equações do Movimento do Corpo Rígido
• Temos finalmente todos os conceitos
especificam o vetor de estados Y(t)
que
 x1 (t ) 


 R (t ) 
Y (t )  
P (t ) 


 L(t ) 
 x(t )   v(t ) 

 


d
d  R(t )   w(t ) R(t ) 
Y (t )  


F (t ) 
dt
dt P(t )

 

 L(t )    (t ) 
46
Codificação Básica em C++
• Assumamos a existência de tipos
denominados de matrix e triple, as quais
implementam operações (soma, subtração
e multiplicação) respectivamente, sobre
matrizes 33 e pontos em 3-d.
47
Estrutura de um Corpo Rígido
struct RigidBody { /* Constant quantities */ double mass; /* mass M */
matrix Ibody, /* Ibody */ Ibodyinv; /* I-1body (inverse of Ibody)
*/
/* State variables */
triple x; /* x(t) */ matrix R; /* R(t) */ triple P, /* P(t) */
L; /* L(t) */
/* Derived quantities (auxiliary variables) */
matrix Iinv; /* I-1(t) */
triple v, /* v(t)*/ omega; /* w(t) */
/* Computed quantities */
triple force, /* F(t)*/ torque; /* τ(t) */
};
/* and assume a global array of bodies */
RigidBody Bodies[NBODIES];
48
Estrutura de um Corpo Rígido
• As quantidades mass, Ibody e Ibodyinv
devem ser previamente calculadas para
cada membro do conjunto de “Bodies”.
• Todas as condições iniciais para cada
corpo rígido também são especificadas
pela atribuição valores às variáveis de
estado x, R, P e L de cada membro de
“Bodies”.
49
Passando parâmetros ao solucionador de
EDO
/* Copy the state information into an array */
void State_to_Array(RigidBody *rb, double *y)
{ *y++ = rb->x[0]; /* x component of position */
*y++ = rb->x[1]; /* etc. */
*y++ = rb->x[2];
for(int i = 0; i < 3; i++) /* copy rotation matrix */
for(int j = 0; j < 3; j++)
*y++ = rb->R[i,j];
*y++ = rb->P[0]; *y++ = rb->P[1];
*y++ = rb->P[2];
*y++ = rb->L[0]; *y++ = rb->L[1];
*y++ = rb->L[2];
}
50
Recebendo parâmetros do solucionador de
EDO
/* Copy information from an array into the state variables */
void Array_to_State(RigidBody *rb, double *y)
{ rb->x[0] = *y++;
rb->x[1] = *y++;
rb->x[2] = *y++;
for(int i = 0; i < 3; i++)
for(int j = 0; j < 3; j++)
rb->R[i,j] = *y++;
rb->P[0] = *y++;
rb->P[1] = *y++;
rb->P[2] = *y++;
rb->L[0] = *y++;
rb->L[1] = *y++;
rb->L[2] = *y++;
rb->v = rb->P / mass;
/* Compute auxiliary variables... */
/* v(t)=P(t)/M */
/* I-1(t)=R(t) I-1body R(t) T*/
}
rb->Iinv = R * Ibodyinv * Transpose(R);
/* w(t)=I-1(t)L(t) */
rb->omega = rb->Iinv * rb->L;
51
Passando e recebendo o estado de todos os
corpos
#define STATE_SIZE 18
void Array_to_Bodies(double y[])
{
for(int i = 0; i < NBODIES; i++)
Array_to_State(&Bodies[i], &y[i * STATE_SIZE]);
}
void Bodies_to_Array(double y[])
{
for(int i = 0; i < NBODIES; i++)
State_to_Array(&Bodies[i], &y[i * STATE_SIZE]);
}
52
Computando dydt
• Suporemos que exista a função que calcula a força F(t) e o Torque τ(t)
agindo sobre cada corpo:
void Compute_Force_and_Torque(double t, RigidBody *rb);
• Assim dydt será:
void dydt(double t, double y[], double ydot[])
{
/* put data in y[] into Bodies[] */
Array_to_Bodies(y);
for(int i = 0; i < NBODIES; i++)
{ Compute_Force_and_Torque(t, &Bodies[i]);
ddt_State_to_Array(&Bodies[i],&ydot[i * STATE_SIZE]);
m}
}
53
Função que Atualiza a Estrutura de dY(t)/dt
do Corpo Rígido
void ddt_State_to_Array(RigidBody *rb, double *ydot)
{
/* copy dx/dt= v(t) into ydot */
*ydot++ = rb->v[0];
*ydot++ = rb->v[1];
*ydot++ = rb->v[2];
/* Compute dR/dt=w(t)* R(t) */
matrix Rdot = Star(rb->omega) * rb->R;
/* copy dR/dt into array */
for(int i = 0; i < 3; i++)
for(int j = 0; j < 3; j++)
*ydot++ = Rdot[i,j];
/*dP/dt=F(t) */
*ydot++ = rb->force[0]; *ydot++ = rb->force[1]; *ydot++ = rb->force[2];
/* dL/dt=τ(t) */
*ydot++ = rb->torque[0];
}
*ydot++ = rb->torque[1]; *ydot++ = rb->torque[2];
54
Rotina Star
• A rotina Star utilizada é definida como:
matrix Star(triple a);
e retorna a seguinte matrix:
 a[2] a[1] 
 0


0
 a[0] 
 a[2]
  a[1] a[0]

0 

55
Executando a Simulação
void RunSimulation()
{ double y0[STATE_SIZE * NBODIES], yfinal[STATE_SIZE * NBODIES];
InitStates();
Bodies_to_Array(yfinal);
for(double t = 0; t < 10.0; t += 1./30.)
{
/* copy yfinal back to y0 */
for(int i = 0; i < STATE_SIZE * NBODIES; i++)
y0[i] = yfinal[i];
ode(y0, yfinal, STATE_SIZE * NBODIES,t, t+1./30., dydt);
/* copy dY(t+1/30.)/t into state variables */
}
}
Array_to_Bodies(yfinal);
DisplayBodies();
56
“Quaternions” x Matrix de Rotações
• A matriz de “quaternions” é um tipo que
comporta apenas 4 elementos;
• Evita-se o uso da matriz de rotação devido a
maior propagação de erros numéricos (“drift”);
• Visualmente
percebe-se
um
efeito
de
“deslizamento”; e
• Enquanto a matriz de rotação faz uso de nove
parâmetros para descrever três graus de
liberdade, os quaternions utilizam quatro
parâmetros com um único descrevendo os três
graus de liberdade
menos “drift” que as
matrizes de rotação.
57
Conclusão
• Foi apresentado, passo a passo, todo o
“background” utilizado na Modelagem
Baseada na Física para implementação de
movimentos de corpos rígidos, assim
como uma pequena codificação em C++
que pode ser utilizada como um “first
step” no desenvolvimento do aplicativo.
58
Bibliografia
• “An Introduction to Physically Based
Modeling: Rigid Body Simulation I –
Unconstrained Rigid Body Dynamics”; David
Baraff (Notas do Curso Siggraph 94/97);
• Física Vol.1; R. Resnick e D. Halliday
• “Physically-Based Modeling for Computer
Graphics”; Ronen Barzel
59