物理引擎如何建模 angular 速度和 angular 加速度（3D）
How do physics engines model angular velocity and angular acceleration (in 3D)
我一直在寻找这个问题的答案，但找不到直接的答案。
好像网上的说法是 angular velocity 存储为一个向量，方向是旋转的轴，长度是旋转的速率。
这有一定道理，但它完全忽略了您实际使用它的方式。
所以有2个关键操作需要完成。
- 应用 angular 加速度（如何将 angular 加速度（大概在相同的数据类型中）添加到您的 angular 速度？）
- 计算刚体上一点的速度，同时考虑angular和线速度。
我研究了像 bullet 这样的物理引擎是如何处理它的，但是数学非常混乱，一些运算符重载似乎隐藏了实际发生的事情。
首选 C++ 或 C 语言的答案，但只要数学清楚，一切都可行。
我建议阅读论文 基于物理的建模简介：刚性 Body 模拟 I - 无约束刚性 Body 动力学 作者 David Baraff。
在论文中，您有刚性 body 的 angular 动量 L 和 angular 速度 ω 的模型。但它不使用 angular 加速度，因为 angular 动量取决于施加于它的扭矩。您还需要定义自己的惯性张量 I(t)，它“描述 body 中的质量如何相对于 body 的质心分布。”在论文中还有一个用 C 语言实现的整个模型（质量、力、扭矩、angular 动量）。
这个 repo 也有一个论文的实现（在 C 中）：RiBOEngine
来源：https://www.cs.cmu.edu/~baraff/sigcourse/notesd1.pdf
在bullet3物理引擎中，实现与上面描述的论文非常相似。如果您查看 btRigidBody::integrateVelocities() method，您会发现为了计算 angular 速度，您需要张量和扭矩来计算它。
数据结构
首先是一些数据结构。
线性代数
我假设有 Vector3 类 用于 3 个值的向量和 Quaternion 类 用于 4 个值的旋转以及 Matrix3 用于 3×3惯性和旋转等操作。每个都遵循其相应的线性代数规则（加法 +、减法 -、缩放 * 和乘积 *）。每个 object 都有自己的规则来处理这些运算符。根据上下文，它可能是比例运算，或 matrix-vector 乘积，或四元数乘积。
Body 属性
每个body都有一个标量质量m，和一个body质量惯性张量I_body 用于运动方程。您还应该保持逆惯性 I_body_inv = inverse( I_body ).
Body 州
body 的当前状态，包含其位置 r (Vector3)，方向 q (Quaternion) 和运动。但通常不是存储两个速度矢量（平移v（Vector3），旋转ω（Vector 3）） , 存储两个动量向量(线性p(Vector 3),angularL(Vector 3 )).每个都在质心处求和。
状态向量有 13 个分量，因为每个向量有 3 个分量，而四元数 q 有 4
Y = [r, q, p, L ]
您可以选择使 Y 成为 array/vector 或 struct 并使用适当的运算符来执行线性代数。
运动提取
给定 Y 你总是可以提取平移速度 v 和 angular 速度 ω[= body 的 134=] 具有以下算法
[v, ω] = GetMotion(Y)
// Get Motion from Y=[r, q, p, L]
r = Y[0:2];
q = Y[3:6];
p = Y[7:9];
L = Y[10:12];
// Get 3×3 rotation matrix R from orientation q
R = rotation(q);
// Get 3×3 inverse mass moment of inertia on the world frame
I_inv = R * I_body_inv * transpose(R);
// Get translational velocity vector
v = p / m;
// Get rotational velocity vector
ω = I_inv * L;
return [v, ω];
还需要反向设置body
的初始条件
Y = GetState(r, q, v, ω)
// Get state vector from position & motion
// Get 3×3 rotation matrix R from orientation q
R = rotation(q);
// Get 3×3 mass moment of inertia on the world frame
I = R * I_body * transpose(R);
// Get translational momentum vector
p = m * v;
// Get angular momentum vector
L = I * ω;
return [r, q, p, L];
整合
在每个时间步，施加的负载向量（力 F、扭矩 τ）都会改变 body 状态向量。
Body 状态率
body状态的时间率，我称Y'，计算如下
Y' = GetRate(t, Y)
// Assemble the time rate of the body state Y at time t
r = Y[0:2];
q = Y[3:6];
p = Y[7:9];
L = Y[10:12];
// Get the motion of the body
[v, ω] = GetMotion(Y);
// Get the loading
[F, τ] = GetLoading(t, r, q, v, ω);
r' = v
q' = 0.5*quaterion(ω)*q
p' = F
L' = τ
return [r', q', p', L'];
其中 quaternion(ω) = [ωx,ωy,ωz,0] 和四元数乘积 * 符合四元数规则。
Runge-Kutta 4
使用 RK4 方法随时间积分很常见
ΔY = GetStep(t, Y, h)
// Integrate state vector Y, over time step h at time t
K0 = GetRate(t, Y);
K1 = GetRate(t+h/2, Y + h/2*K0);
K2 = GetRate(t+h/2, Y + h/2*K1);
K3 = GetRate(t+h, Y + h*K2);
return h/6*(K0 + 2*K1 + 2*K2 + K3);
模拟
模拟是状态向量在许多有限时间步长上的积分。
Body 加载中
需要一个加载函数来描述每个 body 在每个时间帧上的 forces/torques 来评估状态率
[F,τ] = GetLoading(t, r, q, v, ω)
// Describe the force vector F, and torque vector τ
// from time t, position r, orientation q,
// velocity v, and rot. velocity ω
模拟循环
并且模拟循环积分直到达到目标时间
Y = GetState(r, q, v, ω);
time = 0;
while(time<target_time)
{
Y += GetStep(time, Y, h);
time += h;
}
这就是单刚性 body 模拟的基础知识。
Post 处理中
当您拥有 body 状态 Y 时，您将获得 [v, ω] = GetMotion(Y) 的质心运动。但是要回答第二个问题，您需要任意点的速度 v_A。该点位于矢量 r_A
v_A = GetVelocityAtPoint(Y, r_A)
// Get velocity of world point A
r = Y[0:2]
[v, ω] = GetMotion(Y)
return v + cross(ω, r_A - r);
如果你想让一个点v_B的速度骑在body上，并由一个局部向量定义r_B_body 然后
v_B = GetVelocityAtBodyPoint(Y, r_B_body)
// Get velocity of body point B
r = Y[0:2]
q = Y[3:6]
R = rotation(q)
r_A = r + R*r_B_body
return GetVelocityAtPoint(Y, r_A)
请注意，在许多步骤中，方向 Quaternion q 可能会偏离单位四元数，因此需要 re-normlized.
q = normalize(q)
// m = sqrt(qx^2+qy^2+qz^2+qw^2)
// q = q/m
根据我的经验，偏离单位四元数的情况不会逐渐恶化，因为它往往会先偏离然后再回到单位四元数，由此产生的方向角误差实际上非常小。
我一直在寻找这个问题的答案，但找不到直接的答案。 好像网上的说法是 angular velocity 存储为一个向量，方向是旋转的轴，长度是旋转的速率。 这有一定道理，但它完全忽略了您实际使用它的方式。 所以有2个关键操作需要完成。
- 应用 angular 加速度（如何将 angular 加速度（大概在相同的数据类型中）添加到您的 angular 速度？）
- 计算刚体上一点的速度，同时考虑angular和线速度。
我研究了像 bullet 这样的物理引擎是如何处理它的，但是数学非常混乱，一些运算符重载似乎隐藏了实际发生的事情。
首选 C++ 或 C 语言的答案，但只要数学清楚，一切都可行。
我建议阅读论文 基于物理的建模简介：刚性 Body 模拟 I - 无约束刚性 Body 动力学 作者 David Baraff。
在论文中，您有刚性 body 的 angular 动量 L 和 angular 速度 ω 的模型。但它不使用 angular 加速度，因为 angular 动量取决于施加于它的扭矩。您还需要定义自己的惯性张量 I(t)，它“描述 body 中的质量如何相对于 body 的质心分布。”在论文中还有一个用 C 语言实现的整个模型（质量、力、扭矩、angular 动量）。
这个 repo 也有一个论文的实现（在 C 中）：RiBOEngine
来源：https://www.cs.cmu.edu/~baraff/sigcourse/notesd1.pdf
在bullet3物理引擎中，实现与上面描述的论文非常相似。如果您查看 btRigidBody::integrateVelocities() method，您会发现为了计算 angular 速度，您需要张量和扭矩来计算它。
数据结构
首先是一些数据结构。
线性代数
我假设有 Vector3 类 用于 3 个值的向量和 Quaternion 类 用于 4 个值的旋转以及 Matrix3 用于 3×3惯性和旋转等操作。每个都遵循其相应的线性代数规则（加法 +、减法 -、缩放 * 和乘积 *）。每个 object 都有自己的规则来处理这些运算符。根据上下文，它可能是比例运算，或 matrix-vector 乘积，或四元数乘积。
Body 属性
每个body都有一个标量质量m，和一个body质量惯性张量I_body 用于运动方程。您还应该保持逆惯性 I_body_inv = inverse( I_body ).
Body 州
body 的当前状态，包含其位置 r (Vector3)，方向 q (Quaternion) 和运动。但通常不是存储两个速度矢量（平移v（Vector3），旋转ω（Vector 3）） , 存储两个动量向量(线性p(Vector 3),angularL(Vector 3 )).每个都在质心处求和。
状态向量有 13 个分量，因为每个向量有 3 个分量，而四元数 q 有 4
Y = [r, q, p, L ]
您可以选择使 Y 成为 array/vector 或 struct 并使用适当的运算符来执行线性代数。
运动提取
给定 Y 你总是可以提取平移速度 v 和 angular 速度 ω[= body 的 134=] 具有以下算法
[v, ω] = GetMotion(Y)
// Get Motion from Y=[r, q, p, L]
r = Y[0:2];
q = Y[3:6];
p = Y[7:9];
L = Y[10:12];
// Get 3×3 rotation matrix R from orientation q
R = rotation(q);
// Get 3×3 inverse mass moment of inertia on the world frame
I_inv = R * I_body_inv * transpose(R);
// Get translational velocity vector
v = p / m;
// Get rotational velocity vector
ω = I_inv * L;
return [v, ω];
还需要反向设置body
的初始条件Y = GetState(r, q, v, ω)
// Get state vector from position & motion
// Get 3×3 rotation matrix R from orientation q
R = rotation(q);
// Get 3×3 mass moment of inertia on the world frame
I = R * I_body * transpose(R);
// Get translational momentum vector
p = m * v;
// Get angular momentum vector
L = I * ω;
return [r, q, p, L];
整合
在每个时间步，施加的负载向量（力 F、扭矩 τ）都会改变 body 状态向量。
Body 状态率
body状态的时间率，我称Y'，计算如下
Y' = GetRate(t, Y)
// Assemble the time rate of the body state Y at time t
r = Y[0:2];
q = Y[3:6];
p = Y[7:9];
L = Y[10:12];
// Get the motion of the body
[v, ω] = GetMotion(Y);
// Get the loading
[F, τ] = GetLoading(t, r, q, v, ω);
r' = v
q' = 0.5*quaterion(ω)*q
p' = F
L' = τ
return [r', q', p', L'];
其中 quaternion(ω) = [ωx,ωy,ωz,0] 和四元数乘积 * 符合四元数规则。
Runge-Kutta 4
使用 RK4 方法随时间积分很常见
ΔY = GetStep(t, Y, h)
// Integrate state vector Y, over time step h at time t
K0 = GetRate(t, Y);
K1 = GetRate(t+h/2, Y + h/2*K0);
K2 = GetRate(t+h/2, Y + h/2*K1);
K3 = GetRate(t+h, Y + h*K2);
return h/6*(K0 + 2*K1 + 2*K2 + K3);
模拟
模拟是状态向量在许多有限时间步长上的积分。
Body 加载中
需要一个加载函数来描述每个 body 在每个时间帧上的 forces/torques 来评估状态率
[F,τ] = GetLoading(t, r, q, v, ω)
// Describe the force vector F, and torque vector τ
// from time t, position r, orientation q,
// velocity v, and rot. velocity ω
模拟循环
并且模拟循环积分直到达到目标时间
Y = GetState(r, q, v, ω);
time = 0;
while(time<target_time)
{
Y += GetStep(time, Y, h);
time += h;
}
这就是单刚性 body 模拟的基础知识。
Post 处理中
当您拥有 body 状态 Y 时，您将获得 [v, ω] = GetMotion(Y) 的质心运动。但是要回答第二个问题，您需要任意点的速度 v_A。该点位于矢量 r_A
v_A = GetVelocityAtPoint(Y, r_A)
// Get velocity of world point A
r = Y[0:2]
[v, ω] = GetMotion(Y)
return v + cross(ω, r_A - r);
如果你想让一个点v_B的速度骑在body上，并由一个局部向量定义r_B_body 然后
v_B = GetVelocityAtBodyPoint(Y, r_B_body)
// Get velocity of body point B
r = Y[0:2]
q = Y[3:6]
R = rotation(q)
r_A = r + R*r_B_body
return GetVelocityAtPoint(Y, r_A)
请注意，在许多步骤中，方向 Quaternion q 可能会偏离单位四元数，因此需要 re-normlized.
q = normalize(q)
// m = sqrt(qx^2+qy^2+qz^2+qw^2)
// q = q/m
根据我的经验，偏离单位四元数的情况不会逐渐恶化，因为它往往会先偏离然后再回到单位四元数，由此产生的方向角误差实际上非常小。