慣性矩、慣性積、轉動慣量、慣性張量

• 慣性矩是一個幾何量，通常被用作描述截面抵抗彎曲的性質。慣性矩的國際單位為(m⁴)。即面積二次矩，也稱面積慣性矩，而這個概念與質量慣性矩（即轉動慣量）是不同概念。

• 慣性積：質量慣性積是剛體動力學中一個重要的質量幾何性質。剛體中的質量微元 Δm_i與這微元的兩個直角坐標的乘積對剛體的總和。其數值為：$$I_{xy}=\sum_{i}m_ix_iy_i \quad or \quad I_{xy}=\int xydm$$

　　對於給定的物體，慣性積的值與建立的坐標系的位置及方向有關，如果我們選擇的坐標系合適，可使慣性積的值為零。當對於某一坐標軸的慣性積為零時，這種特定的坐標軸稱為慣性主軸或主軸（principal axes），相應的質量慣性矩稱為主慣性矩（principal moments of inertia）。顯然，如果剛體本身具有某種幾何對稱性，那麽它的主軸方向總是沿著它的對稱軸的。但是即使是完全沒有任何對稱性的剛體也是存在慣性主軸的。

• 轉動慣量(Moment of Inertia)是剛體繞軸轉動時慣性的量度。 在經典力學中，轉動慣量又稱質量慣性矩。對於一個質點，J = mr2，其中 m 是其質量，r 是質點和轉軸的垂直距離。

• 慣性張量(Inertia tensor)是描述剛體作定點轉動時的轉動慣性的一組慣性量，其表現形式為由9個分量構成的對稱矩陣，剛體作定點轉動的力學情況要比定軸轉動復雜得多。

　　慣性張量描述了物體的質量分布。以{A}為參考坐標系，剛體相對於坐標系{A}的慣性張量用3×3的矩陣表達：

　　慣性張量矩陣的特征值是主慣性矩，其對應的特征向量是慣性主軸的方向向量（The eigenvalues of an inertia tensor are the principal moments for the body. The associated eigenvectors are the principal axes.） 已知某參考坐標系下的慣性張量，要計算慣性主軸可以參考

牛頓-歐拉方程

　　剛體運動 = 質心的平動 + 繞質心的轉動。其中，質心平動用牛頓方程描述；繞質心的轉動用歐拉方程描述，它們都涉及到質量及其分布。

　　力F作用在剛體的質心上，使得質量為m的剛體產生加速度$\dot{v_c}$，根據牛頓第二運動定律有：$$F=m \dot{v_c}$$

　　作用在剛體上的力矩N使剛體以角速度$\omega$、角加速度$\dot{\omega}$旋轉。根據歐拉方程有：$$N=^{C}I\dot{\omega}+\omega \times ^{C}I \omega$$

　　其中$^{C}I$是以質心參考系{C}為參照的剛體慣性張量。

　　需要註意的是，剛體繞定軸轉動時，角速度矢量$\omega$和角加速度矢量$\dot{\omega}$都是沿著固定的軸線；而剛體繞定點運動時，角速度矢量$\omega$的大小和方向都將不斷變化。角加速度矢量$\dot{\omega}$沿$\omega$的矢端曲線的切線，在一般情況下，它不與角速度矢量$\omega$重合。

剛體定點運動時的動量矩

　　剛體運動中，若剛體或其延拓空間上有一點的位置始終保持不動，這種運動稱為剛體定點運動。剛體定點運動是剛體運動較為復雜的一種運動形式，在工程技術上有著廣泛的應用。定點運動剛體的任何位移都可由繞過定點的某一軸（轉動瞬軸）的一次轉動來實現。瞬軸通過剛體的定點，在不同瞬時有不同的方位。定點運動剛體每瞬時的真實運動就是繞過定點的一系列不同方位的瞬軸的轉動過程。

　　剛體定點運動時，對點O的動量矩為：$$\vec{L}=\sum_{i=1}^{n}\vec{r_i}\times m_i \vec{v_i}$$

　　因為$\vec{v_i}=\vec{\omega} \times \vec{r_i}$

　　故有$$\vec{L}=\sum_{i=1}^{n} m_i[\vec{r_i}\times (\vec{\omega} \times \vec{r_i})]$$

　　上式中各矢量向坐標軸投影，有：

$$\left\{\begin{matrix}
\vec{r_i}=x_i \vec{i}+y_i \vec{j}+z_i \vec{k}\\
\vec{\omega}=\omega_x \vec{i}+\omega_y \vec{j}+\omega_z \vec{k}\\
\vec{L}=L_x \vec{i}+L_y \vec{j}+L_z \vec{k}\
\end{matrix}\right.$$

　　將上面兩式子結合起來寫成矩陣形式如下：

$$\begin{bmatrix}L_x\\ L_y\\ L_z\end{bmatrix}=\begin{bmatrix}
I_{xx} & -I_{xy} & -I_{xz}\\
-I_{yx} & I_{yy} & -I_{yz}\\
-I_{zx}& -I_{zy} & I_{zz}
\end{bmatrix}\begin{bmatrix}
\omega_x\\
\omega_y\\
\omega_z
\end{bmatrix}$$

　　式子中$I_{xx}=\sum m_i(y^2_i + z^2_i)$，$I_{yy}=\sum m_i(x^2_i + z^2_i)$，$I_{zz}=\sum m_i(y^2_i +x^2_i)$，$I_{xy}=I_{yx}=\sum m_i x_i y_i$，$I_{xz}=I_{zx}=\sum m_i x_i z_i$，$I_{yz}=I_{zy}=\sum m_i z_i y_i$

　　簡寫成向量形式為：$$\vec{L}=I \vec{\omega}, \quad I =\begin{bmatrix}
I_{xx} & -I_{xy} & -I_{xz}\\
-I_{yx} & I_{yy} & -I_{yz}\\
-I_{zx}& -I_{zy} & I_{zz}
\end{bmatrix} $$

歐拉動力學方程推導

　　我們首先在靜止坐標系$Oxyz$（固定坐標系/世界坐標系/絕對坐標系）中描述剛體的動力學。剛體有6個自由度，可以選為質心的坐標（3個）和固連在剛體上的動坐標系相對於靜止坐標系的三個歐拉角來表征。假設N為作用在剛體上的總力矩，這裏僅僅需要考慮外力，因為內部的內力的力矩為零。求剛體一般運動的關鍵是應用角動量定力，由已知受力情況去求剛體轉動的角速度：

$$\mathbf{N}= \frac{d\mathbf{L}}{dt}=\frac{I\omega}{dt}$$

　　由於$I$是一個慣性張量，而且要考慮$I\omega$的變化率，如果選固定坐標系作參考，那麽剛體相對固定坐標系的慣性張量$I$的各個分量隨時間改變，這樣求解上式並不容易。我們可以選取剛體的慣量主軸作為參考坐標系（固定在剛體上隨著剛體運動的參照系）使得剛體的慣性張量具有對角形式，且對角線上的分量不隨時間變化。對於這個主軸坐標系，剛體的質量分布是固定的。

　　下面我們試圖在動坐標系（隨體坐標系$O{x}‘{y}‘{z}‘$）中討論剛體的動力學問題。${i}‘$、${j}‘$、${k}‘$為隨動體參考系坐標軸的三個單位矢量，$\omega_{{x}‘}$、$\omega_{{y}‘}$、$\omega_{{z}‘}$為角速度沿三個坐標軸的分量。若$O{x}‘$、$O{y}‘$、$O{z}‘$為三個慣量主軸，則慣性張量中所有慣性積為零。剛體動量矩為：$$L=I_{{x}‘}\omega_{{x}‘}{i}‘+I_{{y}‘}\omega_{{y}‘}{j}‘+I_{{z}‘}\omega_{{z}‘}{k}‘$$

　　可以看出，$L$和$\omega$通常不共線，因為三個轉動慣量一般並不相等，僅當剛體繞慣量主軸轉動時，$L$和$\omega$才共線。

　　下面涉及到兩個不同坐標系觀測到的向量對時間求導的關系 (time derivative in rotating reference frame)。假設剛體以角速度$\omega$繞著某一瞬時軸旋轉，與剛體固連的動坐標系的單位向量$\vec{u}$隨時間變化，變化規律如下：$$\frac{d \vec{u}}{dt}=v=\omega \times \vec{u}$$

　　如果我們在動坐標系中觀測到向量$f$（這時動坐標系的基矢${i}‘$、${j}‘$、${k}‘$相對觀察者是不變的），$f(t)=f_x(t){i}‘+f_y(t){j}‘+f_z(t){k}‘$，現在我們要以靜止參考系為參照（這時${i}‘$、${j}‘$、${k}‘$是隨時間變化的）對$f$求導$$\begin{align*}
\frac{d \mathbf{f}}{dt} &=\frac{df_x}{dt}{i}‘+\frac{d{i}‘}{dt}f_x+ \frac{df_y}{dt}{j}‘+\frac{d{j}‘}{dt}f_y + \frac{df_z}{dt}{k}‘+\frac{d{k}‘}{dt}f_z \\
&=\frac{df_x}{dt}{i}‘+\frac{df_y}{dt}{j}‘+\frac{df_z}{dt}{k}‘+[\mathbf{\omega} \times(f_x {i}‘+f_y {j}‘+f_z {k}‘)]\\
&=(\frac{d \mathbf{f}}{dt})_r+ \mathbf{\omega} \times \mathbf{f}(t)
\end{align*}$$

　　牛頓第一定律和第二定律只適用於慣性參考系，對於非慣性參考系是不適用的。在非慣性參考系中動力學的基本方程不同於慣性系。根據上面的結論，可以推導出隨體坐標系中的歐拉動力學方程：$$\begin{align*}
N=\frac{dL}{dt}&=(\frac{dL}{dt})_r+\omega \times L \\
&=\frac{d(I\omega)}{dt}+\omega \times I\omega\\
&=I\frac{d\omega}{dt}+\omega \times I\omega
\end{align*}$$

　　下標$r$表示向量在剛體上的動坐標系中描述，一般為慣性主軸坐標系。寫成分量方程如下：$$\left\{\begin{matrix}
I_{{x}‘} \frac{d\omega_{{x}‘}}{dt}+(I_{{z}‘}-I_{{y}‘})\omega_{{y}‘}\omega_{{z}‘}=N_{{x}‘}\\
I_{{y}‘} \frac{d\omega_{{y}‘}}{dt}+(I_{{x}‘}-I_{{z}‘})\omega_{{x}‘}\omega_{{z}‘}=N_{{y}‘}\\
I_{{z}‘} \frac{d\omega_{{z}‘}}{dt}+(I_{{y}‘}-I_{{x}‘})\omega_{{x}‘}\omega_{{y}‘}=N_{{z}‘}
\end{matrix}\right.$$

　　其中，$N_{{x}‘}$、$N_{{y}‘}$、$N_{{z}‘}$分別為外力對$O{x}‘$、$O{y}‘$、$O{z}‘$軸之矩，$I_{{x}‘}$、$I_{{y}‘}$、$I_{{z}‘}$為剛體相對過質心的三個慣量主軸的轉動慣量。這就是求解剛體定點運動的基本方程（一般為非線性微分方程組），稱為歐拉方程。需要註意的是，歐拉方程中各矢量的分量都是按照剛體慣性主軸投影的分量。

　　一般來說，不能從上式求出世界坐標系下角速度的三個分量，因為歐拉動力學方程中所取的主軸坐標系是固定在剛體上的。還需要結合歐拉運動學方程（動坐標系中角速度與歐拉角的關系）：$$\left\{\begin{matrix}\begin{align*}
\omega_{{x}‘}&=\dot{\psi}sin\theta sin\phi+\dot{\theta}cos\phi\\
\omega_{{y}‘}&=\dot{\psi}sin\theta cos\phi-\dot{\theta}sin\phi\\
\omega_{{z}‘}&=\dot{\psi}cos\theta +\dot{\phi}
\end{align*}\end{matrix}\right.$$

　　將歐拉動力學方程和運動學方程結合起構成求解定點運動的微分方程組，消去$\omega_{{x}‘}$、$\omega_{{y}‘}$、$\omega_{{z}‘}$後，便得到關於歐拉角$\psi$、$\phi$、$\theta$的二階微分方程。這是一組非線性常微分方程組，它只在特殊條件下才存在解析解。一般情況下，只能通過數值計算得到方程的數值解。

VREP中定義物體質量及質量分布

　　在物體屬性對話框中點擊Show dynamic properties dialog按鈕，彈出剛體動力學屬性對話框。在這個界面中可以定義物體的質量以及質量分布參數：

• Compute mass & inertia properties for the selected shapes: by clicking this button, you can automatically compute the mass and inertia properties of the selected convex shapes, based on a material uniform density. 點擊該按鈕會彈出一個對話框，用戶可以自定義密度大小，接著點擊OK，軟件會根據凸殼的形狀自動計算質量和慣性屬性。比如新建一個默認大小的立方體，體積為0.001m³，密度為500Kg/m³，那麽軟件自動計算物體的質量將會是0.5Kg. 需要註意，不是所有類型的shape都可以自動計算質量屬性，必須是convex類或pure類的shape才行。

• Mass: the mass of the shape. Selected shapes can have their masses easily increased or decreased by a factor 2 with the M=M*2 (for selection) and the M=M/2 (for selection) buttons. This is convenient to quickly find stable simulation parameters by trial-and-error. 除了自動計算外，還可以直接在Mass文本框中輸入物體質量。右邊的兩個按鈕M=M*2 (for selection) 和M=M/2 (for selection)可以很方便的將輸入的質量加倍或減半，調試時很方便。
• Principal moments of inertia / mass: the mass-less (i.e. divided by the mass of the shape) principal moments of inertia. 這個地方可以定義物體繞慣性主軸的（轉動慣量/質量）的值，註意要除以物體質量。

　　通常，規則幾何體繞質心坐標系的轉動慣量可以查表得到，省去了計算過程。在VREP中添加一個默認大小的圓柱體，底面圓的半徑為0.05m，高為0.1m，默認密度為1000Kg/m³. 那麽：

$$\begin{align*}
m &= \rho \pi r^2 h = 0.785398 \\
I_{xx}m^{-1} &= \frac{m}{12} \left(h^2+3 r^2\right)/m=0.00145833 \\
I_{yy}m^{-1} &= \frac{m}{12} \left(h^2+3 r^2\right)/m=0.00145833 \\
I_{zz}m^{-1} &= \frac{m}{2} r^2/m=0.00125
\end{align*}$$

　　下圖分別是長方體、圓柱體、橢球體以質心坐標系為參考的質量慣性矩（轉動慣量）計算公式：

• Pos./orient. of inertia frame & COM relative to shape frame: the configuration of the inertia frame and center of mass expressed relative to the shape‘s reference frame (always located at the geometric center of the shap). 這一選項可以用來更改慣性參考坐標系在shape‘s reference frame中的位置和姿態。shape reference frame是一個原點固定的坐標系，始終位於幾何體中心。在場景中新建一個圓柱體，可以看到紅、綠、藍三個軸組成的坐標系就是shape reference frame，默認情況下inertia frame是與其重合的。我們可以修改一下坐標系的位置（相當於改變了質心位置），將X分量修改為0.05m，即質心偏離中心跑到了圓柱面上：　

　　進行動力學仿真，可以看出圓柱體的質心由於偏離中心太遠，在傾覆力矩的作用下倒下：

• Set inertia matrix & COM relative to absolute frame: 點擊這個按鈕會彈出一個對話框，用戶可以在世界坐標系下（以VREP中的絕對坐標系為參考，而不是部件的局部坐標系）自定義質心位置，並且可以設置剛體的慣性張量/質量矩陣（註意這是一個對稱矩陣，矩陣參數的設置還是以質心坐標系為參考）。

• Inertia matrix divided by the mass: the inertia matrix or tensor. Values are mass-less (i.e. divided by the mass of the shape). The matrix must be expressed relative to the center of mass of the shape (i.e the matrix is symmetric). 設置質心坐標系下的慣性張量/質量矩陣。
• Position of the center of mass: 設置質心在世界坐標系中的位置。
• Apply to selected shapes: when checked, then all selected shapes will have the same inertia properties relative to the absolute reference frame (i.e. all center of masses and inertia matrices will be coincident).

Solidworks中測量質量屬性

　　大多數三維建模軟件（例如 SolidWorks、CATIA）以及一些仿真軟件（如 Adams、VREP）都提供慣性計算功能。下面以SolidWorks為例對介紹一下如何獲取慣性參數，並與VREP進行對比。在Solidworks中新建一個零件，該零件由兩個圓柱體堆疊而成（零件原點在大圓柱體底面中心），半徑分別為50mm、30mm，兩個圓柱體的高都為50mm。為了軟件能正確計算慣性參數，需要給零件設定材料，這裏選為普通碳鋼，密度為7800Kg/m³

　　在工具→評估→質量屬性中可以查看零件的慣性參數（如果單位或顯示的數值不合適可以在選項中調整）

　　註意以下幾點：

• 重心/質心坐標是相對於零件原點的
• SolidWorks列出了3個慣性張量，它們之間的區別就在於分別相對於不同的坐標系：

1. 相對於主軸坐標系；其中的 I_x、I_y、I_z 三個向量表示主軸坐標系相對於繪圖坐標系（零件參考坐標系）的姿態，即主軸坐標系的 x、y、z 三個軸向量在繪圖坐標系中的表示。而 Px、Py、Pz 表示慣性主矩

   $$I_{A}=\left[
   \begin{matrix}
   \rm Px & 0 & 0\\
   0 & \rm Py & 0\\
   0 & 0 & \rm Pz\\
   \end{matrix}
   \right]$$

2. 相對於原點與主軸坐標系重合，但是各軸與繪圖坐標系一致的坐標系
3. 相對於繪圖坐標系（SolidWorks中稱為輸出坐標系）

• 由於該零件有很強的對稱性（以繪圖坐標系為例，零件在Y方向和Z方向對稱），因此慣性張量中除對角線外元素都為0

　　下面將Solidworks中的零件導出成STL網格模型，然後再導入VREP中。這種導入的STL模型在VREP中為Simple random shape：can represent any mesh. It has one color and one set of visual attributes. Not optimised nor recommended for dynamics collision response calculation (since very slow and unstable). 但是由於這種類型的幾何體不支持自動計算質量屬性，因此先將其轉換為凸殼或組合凸殼：

　　然後設置密度為7800Kg/m³，點擊按鈕進行自動計算：

　　經驗算，零件質量、質心位置以及慣性主矩與Solidworks的計算結果一致。區別是VREP中零件的固定參考坐標系是默認在包圍盒的中心位置，目前還沒發現可以更改這一參考點的地方。

參考：

物理引擎中的剛體動力學

基於Mathematica的機器人仿真環境（機械臂篇）

機器人動力學--牛頓-歐拉方程

[常見幾何體]轉動慣量公式表

慣性矩與慣性張量的關系

Cross product_Wikipedia

Introduction to Robotics - Mechanics and Control Chapter 6. Manipulator dynamics

Modern Robotics Mechanics, Planning, and Control Chapter 8. Dynamics of Open Chains

剛體質量分布與牛頓-歐拉方程