刚体姿态动力学推导

找不到服务器zhn / 2023-09-06 / 原文

摘要首先推导角速度公式和角加速度公式，并举一个例子说明角速度公式的细节，然后从加速度公式出发推导转动惯量矩阵 J 和姿态动力学方程 Jw'=w×Jw，并解释为什么需要进行一次坐标变换。
（未完成）

旋转坐标系下的速度和加速度

（角速度公式）旋转坐标系 $\mathcal{F}_2$ 绕惯性坐标系 $\mathcal{F}_1$ 以角速度 $\vec{\omega}$ 旋转，设 $\mathcal{F}_1$ 下位置向量 $\vec r$ 的一阶导和二阶导分别为 $\dot{\vec{r}}_1$ 和 $\ddot{\vec{r}}_1$，$\mathcal{F}_2$ 下位置向量的一阶导和二阶导分别为 $\dot{\vec{r}}_2$ 和 $\ddot{\vec{r}}_2$，$\mathcal{F}_2$，满足

\[\begin{aligned} \dot{\vec{r}}_1 =& \dot{\vec{r}}_2 + \vec{\omega}\times\vec{r}_1 \\ \ddot{\vec{r}}_1 =& \ddot{\vec{r}}_2 + 2\vec{\omega}\times\dot{\vec{r}}_2 + \dot{\vec{\omega}}_2\times\vec{r} + \vec{\omega}\times(\vec{\omega}\times\vec{r}) \end{aligned}\]

证明： 设同一向量 $\vec{r}$ 在坐标系 $\mathcal{F}_1$ 和 $\mathcal{F}_2$ 下的坐标分别为

\[\begin{aligned} \vec{r} =& \begin{bmatrix} \vec{e}_x & \vec{e}_y & \vec{e}_z \end{bmatrix} \begin{bmatrix} x_1 \\ y_1 \\ z_1 \end{bmatrix} = \begin{bmatrix} \vec{a}_x & \vec{a}_y & \vec{a}_z \end{bmatrix} \begin{bmatrix} x_2 \\ y_2 \\ z_2 \end{bmatrix} \\ =& \mathcal{F}_1\vec{r}_1 = \mathcal{F}_2\vec{r}_2 \end{aligned}\]

其中 $[\vec{e}_x\ \vec{e}_y\ \vec{e}_z]$ 表示惯性坐标系 $\mathcal{F}_1$ 下的三轴单位向量，$[\vec{a}_x\ \vec{a}_y\ \vec{a}_z]$ 表示坐标系 $\mathcal{F}_2$ 的三轴单位向量在惯性系 $\mathcal{F}_1$ 下的坐标，
三个单位向量张成旋转坐标系 $\mathcal{F}_2$ ，
则向量 $\vec{r}$ 的一阶导

\[\begin{aligned} \dot{\vec{r}}_1 =& \frac{\text{d}}{\text{d}t}(\mathcal{F}_2\vec{r}_2) \\ =& \mathcal{F}_2\dot{\vec{r}}_2 + \dot{\mathcal{F}_2}\vec{r}_2 \\ =& \mathcal{F}_2\dot{\vec{r}}_2 + \vec{\omega}_1 \times \mathcal{F}_2\vec{r}_2 \\ =& \mathcal{F}_2\dot{\vec{r}}_2 + \vec{\omega}_1 \times \vec{r}_1 \end{aligned}\]

此处以及后面的 $\vec\omega_1$ 都是在惯性系 $\mathcal{F}_1$ 下的坐标，下面简记作 $\vec\omega$。
此处补充一些个人理解。$\mathcal{F}_2=\vec{\omega} \times \mathcal{F}_2$ 是因为 $\mathcal{F}_2$ 可以简单地看作三个向量绕旋转轴 $\vec\omega$ 转动，不涉及坐标变换；而 $\dot{\vec{r}}_1\neq\vec{\omega} \times \vec{r}_1$ 是因为此时向量 $\vec r$ 不是绕旋转轴 $\vec\omega$ 转动，而是个复合运动，$\vec\omega$ 是本体系（旋转坐标系）相对于惯性系的旋转角速度，$\vec r$ 在本体系下还有 $\vec r_2$ 的角速度，当 $\vec r$ 在本体系下的角速度为0时两种情况等价，可以理解成有个本体系固定在 $\vec r$ 上。
在这里插入图片描述
举个例子，如图所示，惯性系用 $Ox_1y_1z_1$ 表示，本体系用 $Ox_2y_2z_2$ 表示。一开始时，本体系与惯性系重合，向量 $\vec r$ 与y轴正向重合，然后本体系沿惯性系x轴逆时针旋转，角速度为 $\vec\omega_1$；向量 $\vec r$ 沿本体系z轴顺时针旋转，角速度为 $\vec\omega_2$

\[\vec\omega_1=\begin{bmatrix} 1 \\ 0 \\ 0 \end{bmatrix},\ \vec\omega_2=\begin{bmatrix} 0 \\ 0 \\ -1 \end{bmatrix}\]

注意两个角速度向量分别在惯性系和本体系下表示。向量 $\vec r$ 在本体系下的坐标为

\[\vec r_2=\begin{bmatrix} \sin t \\ \cos t \\ 0 \end{bmatrix} \]

本体系到惯性系的旋转矩阵为

\[ \mathcal F_2=\begin{bmatrix} 1 & 0 & 0 \\ 0 & \cos t & -\sin t \\ 0 & \sin t & \cos t \end{bmatrix},\ \mathcal F_2^{-1}=\begin{bmatrix} 1 & 0 & 0 \\ 0 & \cos t & \sin t \\ 0 & -\sin t & \cos t \end{bmatrix} \]

惯性系到惯性系的旋转矩阵 $\mathcal F_1$ 为单位阵。所以向量 $\vec r$ 在惯性系下的坐标为

\[ \vec r_1=\mathcal F_2\vec r_2= \begin{bmatrix} 1 & 0 & 0 \\ 0 & \cos t & -\sin t \\ 0 & \sin t & \cos t \end{bmatrix} \begin{bmatrix} \sin t \\ \cos t \\ 0 \end{bmatrix} =\begin{bmatrix} \sin t \\ \cos^2 t \\ \sin t\cos t \end{bmatrix} \]

两边同时求导得到角速度公式

\[\begin{aligned} \dot{\vec{r}}_1 =& \mathcal{F}_2\dot{\vec{r}}_2 + \vec{\omega}_1 \times \vec{r}_1 \\ \begin{bmatrix} \cos t \\ -\sin 2t \\ \cos 2t \end{bmatrix} =& \begin{bmatrix} 1 & 0 & 0 \\ 0 & \cos t & -\sin t \\ 0 & \sin t & \cos t \end{bmatrix} \begin{bmatrix} \cos t \\ -\sin t \\ 0 \end{bmatrix} +\begin{bmatrix} 1 \\ 0 \\ 0 \end{bmatrix} \times\begin{bmatrix} \sin t \\ \cos^2 t \\ \sin t\cos t \end{bmatrix} \\ =& \begin{bmatrix} \cos t \\ -\sin t\cos t \\ -\sin^2t \end{bmatrix} +\begin{bmatrix} 0 \\ -\sin t\cos t \\ \cos^2 t \end{bmatrix} \\ =& \begin{bmatrix} 1 & 0 & 0 \\ 0 & \cos t & -\sin t \\ 0 & \sin t & \cos t \end{bmatrix} \begin{bmatrix} \cos t \\ -\sin t \\ \cos t \end{bmatrix} \end{aligned}\]

对上式进一步求导得

\[\begin{aligned} \ddot{\vec{r}}_1 =& \vec{\omega} \times \mathcal{F}_2\dot{\vec{r}}_2 + \mathcal{F}_2\ddot{\vec{r}}_2 + \dot{\vec{\omega}} \times \mathcal{F}_1\vec{r}_1 \\ &+ \vec{\omega} \times (\vec{\omega} \times \mathcal{F}_2\vec{r}_2 + \mathcal{F}_2\dot{\vec{r}}_2) \\ =& \mathcal{F}_2\ddot{\vec{r}}_2 + 2\vec{\omega} \times \mathcal{F}_2\dot{\vec{r}}_2 \\ &+ \dot{\vec{\omega}} \times \mathcal{F}_1\vec{r}_1 + \vec{\omega} \times (\vec{\omega} \times \mathcal{F}_1\vec{r}_1) \\ =& \ddot{\vec{r}}_2 + 2\vec{\omega}\times\dot{\vec{r}}_2 + \dot{\vec{\omega}}\times\vec{r} + \vec{\omega}\times(\vec{\omega}\times\vec{r}) \end{aligned}\]

姿态动力学方程

\[M=J\dot{\vec w}+\vec w\times J\vec w \]

下面的推导来自文末参考资料[1]。
角动量定理

\[ \frac{\text dH}{\text dt}=M\quad \left(\frac{\text d(mv)}{\text dt}=F\right) \]

其中 $H,M$ 分别为刚体的角动量和作用在刚体上的力矩（括号中与直线运动的质点动量定理对比便于理解）。刚体角动量是刚体内所有质量微元相对基准点 $S$ 的角动量之和

\[ \vec h=m\vec r\times\vec v,\quad H=\iiint\limits_V\vec r\times\dot{\vec r}\text dm \]

其中 $\text dm$ 为刚体的质量微元。

\[\begin{aligned} & \vec r\times\vec v = \vec r\times(\vec w\times\vec r) = r^2\vec w-(\vec r\cdot\vec w)\vec r \\ =& [(x^2+y^2+z^2)w_x-(xw_x+yw_y+zw_z)x]\vec e_x \\ +& [(x^2+y^2+z^2)w_y-(xw_x+yw_y+zw_z)y]\vec e_y \\ +& [(x^2+y^2+z^2)w_z-(xw_x+yw_y+zw_z)z]\vec e_z \\ =& [(y^2+z^2)w_x-xyw_y-xzw_z]\vec e_x \\ +& [-xyw_x+(x^2+z^2)w_y-yzw_z]\vec e_y \\ +& [-xzw_x-yzw_y+(x^2+y^2)w_z]\vec e_z \\ =&\begin{bmatrix} y^2+z^2 & -xy & -xz \\ -xy & x^2+z^2 & -yz \\ -xz & -yz & x^2+y^2 \\ \end{bmatrix} \begin{bmatrix} w_x \\ w_y \\ w_z \end{bmatrix} \\ \end{aligned}\]

于是

\[H=\iiint\limits_V\vec r\times\vec v\text dm=J\vec w \]

其中 $J$ 为刚体的转动惯量矩阵，形式为

\[J=\begin{bmatrix} J_x & -J_{xy} & -J_{xz} \\ -J_{xy} & J_y & -J_{yz} \\ -J_{xz} & -J_{yz} & J_z \\ \end{bmatrix}\]

其中

\[J_x=\iiint\limits_V(y^2+z^2)\text dm \]

其它项类似。
这里需要注意区分惯性系和与刚体固连的本体系，也就是为什么需要进行一次坐标变换。在惯性系中，刚体上各点的坐标不是定值，只有在本体系下的转动惯量矩阵为定值，所以这里需要转换坐标系。角动量在本体系下表示为 $H_s=J_s\vec w$，$J_s$ 为定值。当本体系以角速度 $\vec\omega$ 绕惯性系旋转时，一个向量 $\vec r$ 在两个坐标系下对时间的导数满足

\[\dot{\vec{r}}_1 = \dot{\vec{r}}_2 + \vec{\omega}\times\vec{r}_1 \]

其中 $\vec r_o$ 和 $\vec r_s$ 分别表示向量 $\vec r$ 在惯性系和本体系下的坐标。所以

\[M_o=\dot H_o=\dot H_s+\vec w\times H_s=J_s\dot{\vec w}+\vec w\times J_s\vec w \]

\[\begin{aligned} & \frac{\text d}{\text dt}(\vec r\times\vec v) \\ =& \dot{\vec r}\times\vec v +\vec r\times\dot{\vec v} \\ =& \vec v\times\vec v +\vec r\times\frac{\text d}{\text dt}(\vec\omega\times\vec r) \\ =& \vec r\times(\dot{\vec\omega}\times\vec r+\vec\omega\times\dot{\vec r}) \\ =& (\vec r\cdot\vec r)\dot{\vec\omega}-(\vec r\cdot\dot{\vec\omega})\vec r +\vec r\times(\vec\omega\times(\omega\times\vec r)) \\ =& r^2\dot{\vec\omega}-(\vec r\cdot\dot{\vec\omega})\vec r +\vec r\times((\vec\omega\cdot\vec r)\vec\omega -(\vec\omega\cdot\vec\omega)\vec r) \\ =& r^2\dot{\vec\omega}-(\vec r\cdot\dot{\vec\omega})\vec r +(\vec\omega\cdot\vec r)(\vec r\times\vec\omega) \end{aligned}\]