刚体状态
刚体的运动状态有两种：平移和旋转。
将某个物体从局部坐标系变化到全局坐标系
p
=
x
+
R
p
0
\mathbf{p} =\mathbf{x} +R\mathbf{p_0}
p=x+Rp0
-
p
0
\mathbf{p_0}
p0 ：物体在局部坐标中的某个点坐标。
- R ：旋转矩阵。
-
x
\mathbf{x}
x ：物体质心在全局坐标系中的坐标（物体在局部坐标中的质心为原点）。
-
p
\mathbf{p}
p ：物体该点在全局坐标中的坐标。
但若我们想知道这个点在世界坐标系中是如何随时间变化的（在这里可以理解为某物体的局部坐标系如何随时间变化到世界坐标系），就需要对时间求导。
对时间求导
p
′
=
x
′
+
R
′
p
0
\mathbf{p}' = \mathbf{x}' + {R} '\mathbf{p_0}
p′=x′+R′p0
因为物体为刚体，所以
p
0
\mathbf {p_0}
p0为定值。可归纳为：
p
′
=
v
+
R
′
p
0
\mathbf{p}' = \mathbf{v} + {R} '\mathbf{p_0}
p′=v+R′p0
对矩阵求导
由于局部坐标系变世界坐标系的R为
R
=
[
u
x
^
,
u
y
^
,
u
z
^
]
R = \left [ \hat{u_x} , \hat{u_y},\hat{u_z} \right ]
R=[ux^,uy^,uz^]
u
x
u_x
ux为局部坐标系x轴在世界坐标系下的坐标。
u
y
,
u
z
u_y,u_z
uy,uz同理。
根据公式：
r
′
=
ω
×
r
{\mathbf{r} }' = \omega \times \mathbf{r}
r′=ω×r
-
r
\mathbf{r}
r 为旋转中心到p的向量
即
R
′
=
[
u
x
^
′
,
u
y
^
′
,
u
z
^
′
]
{R}' = \left [ {\hat{u_x}}' , {\hat{u_y}}' ,{\hat{u_z}}' \right ]
R′=[ux^′,uy^′,uz^′]
为
R
′
=
[
ω
×
u
x
^
,
ω
×
u
y
^
,
ω
×
u
z
^
]
{R}' = \left [ \omega \times \hat{u_x} , \omega \times \hat{u_y} , \omega \times \hat{u_z} \right ]
R′=[ω×ux^,ω×uy^,ω×uz^]
-
ω
为
角
速
度
，
是
个
向
量
\omega 为角速度，是个向量
ω为角速度，是个向量
表示为更简单的形式
R
′
=
ω
∗
R
{R}' = \omega^*R
R′=ω∗R
其中
ω
∗
=
[
0
−
ω
z
ω
y
ω
z
0
−
ω
x
−
ω
y
ω
x
0
]
\omega^* = \begin{bmatrix} 0 & -\omega_z & \omega_y\\ \omega_z & 0 & -\omega_x\\ -\omega_y& \omega_x & 0 \end{bmatrix}
ω∗=⎣⎡0ωz−ωy−ωz0ωxωy−ωx0⎦⎤
惯性
在线性运动中，线动量为
P
\mathbf P
P
P
=
m
v
\mathbf P = m\mathbf v
P=mv
而在旋转运动中，角动量
L
\mathbf L
L
L
=
I
ω
\mathbf L = I \mathbf{\omega}
L=Iω
刚体属性
1.质心
质心的离散计算公式：
c
=
∑
m
i
x
i
∑
m
i
c = \frac{\sum m_i\mathbf{x} _i}{\sum m_i}
c=∑mi∑mixi
计算方法
体素法
将物体分成小体积块，将这些小体积块当作质点。
c
=
1
n
∑
i
=
1
n
x
i
c = \frac{1}{n}\sum_{i=1}^{n} \mathbf{x_i}
c=n1i=1∑nxi
直接计算法
将质心和每个表面三角形连接为四面体。四面体的质心为质点，四面体的体积用来确定质量，最后将各质点按质量加权平均得到总体的质心。
四面体体积
V
=
1
6
d
⋅
(
e
×
f
)
V = \frac{1}{6} \mathbf{d} \cdot (\mathbf{e}\times \mathbf{f} )
V=61d⋅(e×f)
-
d
,
e
,
f
\mathbf d,\mathbf e,\mathbf f
d,e,f为表面三角形三点坐标。
- 总体的质心为（0，0，0）点
四面体的中心
x
=
1
4
(
d
⋅
e
×
f
)
\mathbf{x} = \frac{1}{4} (\mathbf{d} \cdot \mathbf{e} \times \mathbf{f} )
x=41(d⋅e×f)
2.惯性张量
I
=
[
∑
(
y
i
2
+
z
i
2
)
m
i
−
∑
(
x
i
y
i
)
m
i
−
∑
(
x
i
z
i
)
m
i
−
∑
(
x
i
y
i
)
m
i
∑
(
x
i
2
+
z
i
2
)
m
i
−
∑
(
y
i
z
i
)
m
i
−
∑
(
x
i
z
i
)
m
i
−
∑
(
y
i
z
i
)
m
i
∑
(
x
i
2
+
y
i
2
)
m
i
]
I=\left[\begin{array}{lll} \sum\left(y_{i}^{2}+z_{i}^{2}\right) m_{i} & -\sum\left(x_{i} y_{i}\right) m_{i} & -\sum\left(x_{i} z_{i}\right) m_{i} \\ -\sum\left(x_{i} y_{i}\right) m_{i} & \sum\left(x_{i}^{2}+z_{i}^{2}\right) m_{i} & -\sum\left(y_{i} z_{i}\right) m_{i} \\ -\sum\left(x_{i} z_{i}\right) m_{i} & -\sum\left(y_{i} z_{i}\right) m_{i} & \sum\left(x_{i}^{2}+y_{i}^{2}\right) m_{i} \end{array}\right]
I=⎣⎡∑(yi2+zi2)mi−∑(xiyi)mi−∑(xizi)mi−∑(xiyi)mi∑(xi2+zi2)mi−∑(yizi)mi−∑(xizi)mi−∑(yizi)mi∑(xi2+yi2)mi⎦⎤
-
x
i
,
y
i
,
z
i
为
第
i
个
点
的
坐
标
，
m
i
为
第
i
个
点
的
质
量
x_i,y_i,z_i为第i个点的坐标，m_i为第i个点的质量
xi,yi,zi为第i个点的坐标，mi为第i个点的质量
同样，我们可以用多种方法来计算这个矩阵。如体素法，直接计算法等。
世界坐标系中的惯性变量
已知
- 局部坐标系中的惯性变量为
I
0
I_0
I0，为一个固定值。
- R为局部坐标系变为全局坐标系的旋转矩阵。
故世界坐标系中的惯性张量为：
I
=
R
I
0
R
T
I = RI_0R^T
I=RI0RT
因为R为正交矩阵，所以
R
T
=
R
−
1
R^T = R^{-1}
RT=R−1
L
=
I
ω
=
R
I
0
R
T
ω
L = I \omega = RI_0R^T\omega
L=Iω=RI0RTω 可以理解为
- 使用
R
−
1
R^{-1}
R−1将角速度从世界坐标中变化到局部坐标中
- 在局部坐标中使用
I
0
I_0
I0惯性张量计算局部坐标下的角动量。
- 再将局部坐标系下的角动量通过
R
R
R再变回世界坐标系下。
即为
L
=
I
ω
=
R
(
I
0
(
R
T
ω
)
)
L = I \omega = R(I_0(R^T\omega))
L=Iω=R(I0(RTω))
注：在图形学中，可通过将矩阵运算从右向左进行来理解步骤。
刚体运动
力矩
力矩
τ
\tau
τ 度量了角动量
L
L
L 的变化率，有:
τ
=
L
′
\tau = {L}'
τ=L′
力矩的计算与 力臂
r
\mathbf r
r（物体质心到受理点的矢量）和 力
F
\mathbf F
F 有关，有：
τ
=
r
×
F
\tau = \mathbf r \times \mathbf F
τ=r×F
刚体的固定属性
- 质量
m
m
m
- 惯性张量
I
I
I
当前刚体状态
当前刚体状态有如下几个量：
- 线性运动位置
x
\mathbf x
x，为该物体质心从局部坐标运动到世界坐标需要的平移向量
- 角运动
q
或
R
\mathbf q 或R
q或R（
q
\mathbf q
q为四元数），为该物体局部坐标系旋转到世界坐标系下需要的旋转向量或矩阵
- 线动量
P
\mathbf P
P
- 角动量
L
\mathbf L
L
状态导数
状态导数表示刚体状态的变化趋势
- 线性运动的变化趋势为速度
v
v
v
- 角运动的变化趋势
q
˙
\dot{\mathbf q}
q˙ 或
R
˙
\dot{R}
R˙ 为
1
2
ω
q
\frac{1}{2} \omega \mathbf q
21ωq 或
ω
∗
R
\omega^*R
ω∗R.
- 线动量的变化趋势为力
F
F
F
- 角动量的变化趋势为力矩
τ
\tau
τ
相关证明可搜索相关书籍 《基于物理的建模与动画》 或其他资料，这里就不赘述。
通过刚体的固有属性、当前状态以及受力情况计算状态导数
- 线速度
v
v
v
v
=
1
m
P
v = \frac{1}{m}\mathbf P
v=m1P
- 四元数变化率
q
˙
\dot{\mathbf{q}}
q˙ (这里角运动的变化趋势我们用四元数计算)
R
=
R
(
q
)
−
−
−
−
−
−
−
−
−
−
I
−
1
=
R
I
0
−
1
R
T
−
−
−
−
−
−
−
−
−
−
ω
=
I
−
1
L
−
−
−
−
−
−
−
−
−
−
q
˙
=
1
2
ω
q
(
这
里
的
ω
为
四
元
数
(
0
,
ω
)
)
R = R(\mathbf q) \\ ---------- \\ I^{-1} = RI_0^{-1}R^T \\----------\\ \omega = I^{-1}L \\----------\\ \dot{\mathbf{q}} = \frac{1}{2}\omega \mathbf q \\(这里的\omega为四元数(0,\omega))
R=R(q)−−−−−−−−−−I−1=RI0−1RT−−−−−−−−−−ω=I−1L−−−−−−−−−−q˙=21ωq(这里的ω为四元数(0,ω))
- 合力
F
\mathbf F
F
这里需要根据物体所处的环境计算，例如风力，摩擦力，以及之后要说的刚体碰撞等。
- 力矩
τ
\tau
τ
对于每一个非彻体力（只作用在物体的某个点，而不是整个物体的力），有
r
i
=
p
i
−
x
\mathbf {r_i} = \mathbf{p_i} - \mathbf x
ri=pi−x
-
x为物体质心在世界坐标系下的坐标。
-
p
i
\mathbf{p_i}
pi为世界坐标系下受力点的坐标
-
r
\mathbf r
r为力臂。
τ
=
∑
r
i
×
F
i
\tau= \sum \mathbf {r_i} \times \mathbf {F_i}
τ=∑ri×Fi
自此，我们得到了某一时刻状态到下一时刻状态的变化率。现在便可以通过这些变量计算下一时刻状态。
计算下一时刻状态
x
n
e
w
=
x
+
v
△
t
\mathbf x_{new} = \mathbf x + \mathbf v\bigtriangleup{t}
xnew=x+v△t
q
n
e
w
=
q
+
q
˙
△
t
\mathbf q_{new} = \mathbf q + \dot{\mathbf q}\bigtriangleup{t}
qnew=q+q˙△t
P
n
e
w
=
P
+
F
△
t
\mathbf P_{new} = \mathbf P + \mathbf F\bigtriangleup{t}
Pnew=P+F△t
L
n
e
w
=
L
+
τ
△
t
\mathbf L_{new} = \mathbf L + \mathbf \tau\bigtriangleup{t}
Lnew=L+τ△t
故，若将
x
,
q
,
P
,
L
\mathbf x ,\mathbf q,\mathbf P,\mathbf L
x,q,P,L设为一个状态
S
t
a
t
e
State
State
S
S
S，有：
S
n
e
w
=
S
+
S
˙
h
S_{new} = S + \dot{S}h
Snew=S+S˙h
这里
h
h
h表示为步时。
