剛體運動力學02 剛體旋轉運動方程式(Rotational Equations of Motion)

前言

在上一篇剛體運動力學01 中，我們定義了線動量 $\mathbf{p}$、角動量 $\mathbf{h}$ 以及慣量矩陣 $\mathbf{I}$。本篇將進入核心部分：推導剛體的旋轉運動方程式。

這部分的推導之所以複雜，是因為剛體在旋轉時，其慣量矩陣在慣性坐標系下是 隨時間變化 的（$\dot{\mathbf{I}} \neq \mathbf{0}$）。我們將一步步證明這個微分關係，並最終導出適用於飛彈與飛行載具的姿態控制模型。

牛頓第二運動定律 (旋轉部分)

根據牛頓第二運動定律，作用在剛體質心上的力矩總和，等於其角動量的時變率。 若 $\mathbf{m}_i$ 為作用在質心上的合力矩，則有：

\[\sum \mathbf{m}_i = \mathbf{\dot h}\]

將角動量 $\mathbf{h} = \mathbf{I}_c \mathbf{\omega}_c$ 代入（其中 $\mathbf{I}_c$ 為質心慣量矩陣，$\mathbf{\omega}_c$ 為質心角速度）：

\[\begin{aligned} \sum \mathbf{m}_i &= \frac{d}{dt}(\mathbf{I}_c \mathbf{\omega}_c) \\ &= \mathbf{\dot I}_c \mathbf{\omega}_c + \mathbf{I}_c \mathbf{\dot \omega}_c \\ &= \mathbf{\dot I}_c \mathbf{\omega}_c + \mathbf{I}_c \mathbf{\alpha}_c \end{aligned} \tag{1}\]

其中 $\mathbf{\alpha}_c$ 為角加速度。 關鍵問題：由於剛體在空間中旋轉，從世界坐標系（慣性坐標）觀察，剛體的質量分布隨時在改變，因此 $\mathbf{\dot I}_c \neq \mathbf{0}$。接下來我們必須推導 $\mathbf{\dot I}_c$ 的數學形式。

慣量矩陣的微分推導

為了計算 \(\mathbf{\dot I}_c\)，我們引入跟隨剛體旋轉的 體軸坐標系 ($R_{br}$)。 在 $R_{br}$ 中，剛體是靜止的，因此其慣量矩陣 \(\mathbf{I}_{br,c}\) 是常數，即：

\[\mathbf{\dot I}_{br,c} = \mathbf{0}\]

利用旋轉矩陣 \(\mathbf{R}_{br}^0\)（將 $R_{br}$ 轉換至 $R_0$ 的矩陣），我們可以聯繫兩者：

\[\mathbf{I}_c = (\mathbf{R}_{br}^0)^T \mathbf{I}_{br,c} \mathbf{R}_{br}^0\]

註：根據論文符號定義，此處 \(\mathbf{R}_{br}^0\) 為將世界坐標轉至體軸坐標的矩陣，故其轉置 $(\mathbf{R}_{br}^0)^T$ 才是將體軸轉至世界坐標的轉換矩陣。

對上式進行時間微分（利用 Chain Rule）：

\[\mathbf{\dot I}_c = \dot{\mathbf{R}}_{br}^{0^T} \mathbf{I}_{br,c} \mathbf{R}_{br}^0 + \mathbf{R}_{br}^{0^T} \mathbf{I}_{br,c} \dot{\mathbf{R}}_{br}^0 \tag{2}\]

推導旋轉矩陣的導數 $\dot{\mathbf{R}}_{br}^0$： 設矩陣 \(\mathbf{R}_{br}^{0^T}\) 的各行向量為 \((\mathbf{e}_{br,x}, \mathbf{e}_{br,y}, \mathbf{e}_{br,z})\)，代表體軸各軸在世界坐標中的單位向量。其微分等於角速度 $\mathbf{\omega}_c$ 與該向量的外積：

\[\begin{aligned} \dot{\mathbf{R}}_{br}^{0^T} &= \frac{d}{dt} \begin{bmatrix} \mathbf{e}_{br,x} & \mathbf{e}_{br,y} & \mathbf{e}_{br,z} \end{bmatrix} \\ &= \begin{bmatrix} \mathbf{\omega}_c \times \mathbf{e}_{br,x} & \mathbf{\omega}_c \times \mathbf{e}_{br,y} & \mathbf{\omega}_c \times \mathbf{e}_{br,z} \end{bmatrix} \\ &= [\mathbf{\omega}_c \times] \begin{bmatrix} \mathbf{e}_{br,x} & \mathbf{e}_{br,y} & \mathbf{e}_{br,z} \end{bmatrix} \\ &= [\mathbf{\omega}_c \times] \mathbf{R}_{br}^{0^T} \end{aligned}\]

其中 \([\mathbf{\omega}_c \times]\) 為角速度的反對稱矩陣。 同理，對上式取轉置可得 \(\dot{\mathbf{R}}_{br}^0\)：

\[\begin{aligned} \dot{\mathbf{R}}_{br}^0 &= (\dot{\mathbf{R}}_{br}^{0^T})^T \\ &= (\mathbf{R}_{br}^{0^T})^T [\mathbf{\omega}_c \times]^T \\ &= \mathbf{R}_{br}^0 (-[\mathbf{\omega}_c \times]) \\ &= -\mathbf{R}_{br}^0 [\mathbf{\omega}_c \times] \end{aligned}\]

代回 $\mathbf{\dot I}_c$： 將上述導數關係代回式 (2)：

\[\begin{aligned} \mathbf{\dot I}_c &= ([\mathbf{\omega}_c \times] \mathbf{R}_{br}^{0^T}) \mathbf{I}_{br,c} \mathbf{R}_{br}^0 + \mathbf{R}_{br}^{0^T} \mathbf{I}_{br,c} (-\mathbf{R}_{br}^0 [\mathbf{\omega}_c \times]) \\ &= [\mathbf{\omega}_c \times] (\mathbf{R}_{br}^{0^T} \mathbf{I}_{br,c} \mathbf{R}_{br}^0) - (\mathbf{R}_{br}^{0^T} \mathbf{I}_{br,c} \mathbf{R}_{br}^0) [\mathbf{\omega}_c \times] \\ &= [\mathbf{\omega}_c \times] \mathbf{I}_c - \mathbf{I}_c [\mathbf{\omega}_c \times] \end{aligned} \tag{3}\]

這就是慣量矩陣在旋轉坐標系下的微分關係式。

尤拉運動方程式 (Euler Equation)

將推導出的 $\mathbf{\dot I}_c$ (3) 代回最初的牛頓定律 (1)：

\[\begin{aligned} \sum \mathbf{m}_i &= ([\mathbf{\omega}_c \times] \mathbf{I}_c - \mathbf{I}_c [\mathbf{\omega}_c \times]) \mathbf{\omega}_c + \mathbf{I}_c \mathbf{\alpha}_c \\ &= [\mathbf{\omega}_c \times] \mathbf{I}_c \mathbf{\omega}_c - \mathbf{I}_c \underbrace{[\mathbf{\omega}_c \times] \mathbf{\omega}_c}_{\mathbf{\omega} \times \mathbf{\omega} = \mathbf{0}} + \mathbf{I}_c \mathbf{\alpha}_c \\ &= \mathbf{\omega}_c \times (\mathbf{I}_c \mathbf{\omega}_c) + \mathbf{I}_c \mathbf{\alpha}_c \end{aligned}\]

利用 $\mathbf{h} = \mathbf{I}_c \mathbf{\omega}_c$，可得 通用尤拉運動方程式 (General Euler Equation)：

\[\sum \mathbf{m}_i = \mathbf{\omega}_c \times \mathbf{h} + \mathbf{I}_c \mathbf{\alpha}_c \tag{4}\]

此方程式是在 慣性坐標系 (World Frame) 下描述的。但在實際應用（如飛彈控制）中，感測器與致動器都是安裝在機體上的，因此我們需要將其轉換到 體軸坐標系 (Body Frame)。

在等式兩邊同乘旋轉矩陣 $\mathbf{R}_0^{br}$（即 \(\mathbf{R}_{br}^0\)）：

\[\mathbf{R}_0^{br} \sum \mathbf{m}_i = \mathbf{R}_0^{br} (\mathbf{\omega}_c \times \mathbf{I}_c \mathbf{\omega}_c) + \mathbf{R}_0^{br} \mathbf{I}_c \mathbf{\alpha}_c\]

利用向量旋轉性質，可導出體軸下的形式：

\[\sum \mathbf{m}_{br,i} = \mathbf{\omega}_{c,br} \times \mathbf{I}_{c,br} \mathbf{\omega}_{c,br} + \mathbf{I}_{c,br} \mathbf{\alpha}_{c,br} \tag{5}\]

其中下標 $br$ 代表 Body Reference。

主軸簡化與狀態方程式

若選取的體軸坐標系恰好為剛體的 主軸 (Principal Axes)，則慣量矩陣 $\mathbf{I}_{c,br}$ 為對角矩陣： \(\mathbf{I}_{c,br} = \begin{bmatrix} I_x & 0 & 0 \\ 0 & I_y & 0 \\ 0 & 0 & I_z \end{bmatrix}\)

將式 (5) 展開為三個純量方程式：

\[\begin{aligned} M_x &= I_x \dot{\omega}_x + (I_z - I_y) \omega_y \omega_z \\ M_y &= I_y \dot{\omega}_y + (I_x - I_z) \omega_x \omega_z \\ M_z &= I_z \dot{\omega}_z + (I_y - I_x) \omega_x \omega_y \end{aligned} \tag{6}\]

這就是著名的 尤拉運動方程式。

控制系統狀態方程式： 為了設計控制器，我們常將角加速度 $\dot{\omega}$ 移至等號左邊，寫成狀態方程式形式：

\[\begin{bmatrix} \dot{\omega}_x \\ \dot{\omega}_y \\ \dot{\omega}_z \end{bmatrix} = \begin{bmatrix} \frac{(I_y - I_z)\omega_y \omega_z}{I_x} \\ \frac{(I_z - I_x)\omega_z \omega_x}{I_y} \\ \frac{(I_x - I_y)\omega_x \omega_y}{I_z} \end{bmatrix} + \begin{bmatrix} I_x & 0 & 0 \\ 0 & I_y & 0 \\ 0 & 0 & I_z \end{bmatrix}^{-1} \begin{bmatrix} M_x \\ M_y \\ M_z \end{bmatrix}\]

寫成向量形式即為論文中飛彈模型的旋轉動力學部分：

\[\mathbf{\dot \omega} = \mathbf{f}_1(\mathbf{\omega}) + \mathbf{J}^{-1} \mathbf{m}\]

其中 $\mathbf{f}_1(\mathbf{\omega})$ 為科氏力與慣量耦合造成的非線性項，$\mathbf{m}$ 為氣動力矩與控制力矩。

參考資料

1. 林昱廷，飛彈姿態運動之適應控制研究，碩士論文，機械工程學系，國立臺灣科技大學，2021。
2. R. C. Hibbeler and K. B. Yap, Engineering Mechanics: Dynamics, 14th ed., 2017.