剛體質量分佈與牛頓-歐拉方程
慣性矩、慣性積、轉動慣量、慣性張量
- 慣性矩是一個幾何量，一般被用做描述截面抵抗彎曲的性質。慣性矩的國際單位爲(m4)。即面積二次矩，也稱面積慣性矩，而這個概念與質量慣性矩（即轉動慣量）是不一樣概念。
- 慣性積：質量慣性積是剛體動力學中一個重要的質量幾何性質。剛體中的質量微元 Δmi與這微元的兩個直角座標的乘積對剛體的總和。其數值爲：$$I_{xy}=\sum_{i}m_ix_iy_i \quad or \quad I_{xy}=\int xydm$$
對於給定的物體，慣性積的值與創建的座標系的位置及方向有關，若是咱們選擇的座標系合適，可以使慣性積的值爲零。當對於某一座標軸的慣性積爲零時，這種特定的座標軸稱爲慣性主軸或主軸（principal axes），相應的質量慣性矩稱爲主慣性矩（principal moments of inertia）。顯然，若是剛體自己具備某種幾何對稱性，那麼它的主軸方向老是沿着它的對稱軸的。可是即便是徹底沒有任何對稱性的剛體也是存在慣性主軸的。html
- 轉動慣量(Moment of Inertia)是剛體繞軸轉動時慣性的量度。 在經典力學中，轉動慣量又稱質量慣性矩。對於一個質點，J = mr²，其中 m 是其質量，r 是質點和轉軸的垂直距離。
- 慣性張量(Inertia tensor)是描述剛體做定點轉動時的轉動慣性的一組慣性量，其表現形式爲由9個份量構成的對稱矩陣，剛體做定點轉動的力學狀況要比定軸轉動複雜得多。
慣性張量描述了物體的質量分佈。以{A}爲參考座標系，剛體相對於座標系{A}的慣性張量用3×3的矩陣表達：express
慣性張量矩陣的特徵值是主慣性矩，其對應的特徵向量是慣性主軸的方向向量（The eigenvalues of an inertia tensor are the principal moments for the body. The associated eigenvectors are the principal axes.） 已知某參考座標系下的慣性張量，要計算慣性主軸能夠參考慣性主軸計算。less
同一貫量在不一樣參考系中的描述
考慮2個與剛體固連的座標系$B$和$B'$，其中座標系$B$的基向量爲${\vec{b_x},\vec{b_y},\vec{b_z}}$，座標系$B'$的基向量爲${\vec{b_x'},\vec{b_y'},\vec{b_z'}}$。角動量及角速度在這兩個座標系中的描述爲（The angular momentum vector and angular velocity vector of a rigid body are expressed in terms of these basis vectors, as follows ）：dom
能夠看出在不一樣座標系下，描述同一個矢量的座標是不同的（由於參考基準不同）。注意當咱們說起「速度」、「角速度」、「動量矩」等詞時若是沒有指定參考系，則默認這些矢量是在靜止座標系（慣性系）中測量的。物體運動的描述是須要參照物的，好比說剛體角速度爲$\omega$，意味着該剛體相對於靜止參考系的角速度爲$\omega$。ide
角動量和角速度的關係以下：$$\mathbf{L}=\mathbf{I} \mathbf{\omega} \\ \mathbf{L'}=\mathbf{I'} \mathbf{\omega'}$$工具
其中$\mathbf{L}=(L_x,L_y,L_z)$，$\mathbf{L'}=(L'_x,L'_y,L'_z)$，$\mathbf{\omega}=(\omega_x,\omega_y,\omega_z)$，$\mathbf{\omega'}=(\omega'_x,\omega'_y,\omega'_z)$，$\mathbf{I}$是剛體關於座標系$B$的慣性張量，$\mathbf{I'}$是剛體關於座標系$B'$的慣性張量。ui
考慮座標系$B$到$B'$的轉換矩陣爲$R$，則有下列關係成立：$$L'=RL,\omega'=R \omega$$this
結合上面兩個公式，能夠獲得：$$ I'=RIR^{T}$$spa
牛頓-歐拉方程
剛體運動 = 質心的平動 + 繞質心的轉動。其中，質心平動用牛頓方程描述；繞質心的轉動用歐拉方程描述，它們都涉及到質量及其分佈。.net
力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$：$$\Omega=\begin{bmatrix}
0& - \omega_3& \omega_2\\
\omega_3& 0& - \omega_1\\
- \omega_2& \omega_1& 0
\end{bmatrix}$$
利用等式：$\omega \times I \omega=\Omega I \omega$，歐拉方程也能夠簡寫爲：$$I\dot{\omega}+\Omega I\omega=N$$
須要注意的是，剛體繞定軸轉動時，角速度矢量$\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{d(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$才共線。注意公式中$L$、$I$、$\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*}$$
根據上面的結論，能夠推導出歐拉動力學方程。Euler方程就是把動量矩定理表示在動座標系/非慣性系（一般是與剛體固連的座標系中）：$$\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}'}$爲剛體相對過質心的三個慣量主軸的轉動慣量。這就是求解剛體定點運動的基本方程（通常爲非線性微分方程組），稱爲歐拉方程。上面歐拉方程中各矢量的份量都是按照剛體慣性主軸投影的份量。公式中的物理量與參考系的關係很容易弄混淆，能夠參考Rigid-Body Dynamics中關於歐拉方程的詳細描述，也能夠參考下面幾個問題：
1. Euler equation and conservation of angular momentum (rigid body) ：$\omega$ in Euler's equations refer to the angular velocity vector expressed in the (moving) body frame of reference. And because the frame of reference is moving the description of the vector is non-constant.
2. Why is body frame angular velocity nonzero? ： Angular velocity is being measured with respect to an inertial frame, but its components can be taken with respect to any basis we wish such as one rotating with the body.
關於物理量與參考座標系之間關係的描述也能夠參考機器人經典書籍：Modern Robotics Mechanics, Planning, and Control. Chapter 3 Rigid-Body Motions. A Word about Vectors and Reference Frames
若是施加在剛體上的力矩是Space-Fxied的，即相對全局座標系固定，不隨剛體運動而變化，那麼就不能將其直接帶入上面歐拉公式的右端，須要轉換到局部參考座標系下才能夠。若是是Body-Moving的力矩，即做用在剛體的局部座標系中，隨着剛體運動，那麼就能夠直接帶入上面的歐拉公式。
須要注意的是：上面Euler方程中各量均爲在與剛體固連的座標系中表達，但這並不意味着此方程只在與剛體固連的座標系中才成立。事實上很容易證實：對任一座標系均有$$I \dot{\omega}+\omega \times I \omega =N$$
式中各量均爲在上述給定座標系中的表達式。固然，將Euler方程表示在與剛體固連座標系中尤爲特殊的優勢，即慣性張量矩陣是常值矩陣，它能夠預先計算或辨識出來，這給Euler方程的應用帶來很大方便。若是歐拉方程在慣性座標系中描述，那麼慣性張量是個隨着剛體運動而變化的變量，須要實時計算該值。
通常來講，不能從上式求出世界座標系下角速度的三個份量，由於歐拉動力學方程中所取的主軸座標系是固定在剛體上的。還須要結合歐拉運動學方程（動座標系中角速度與歐拉角的關係）：$$\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.001m3，密度爲500Kg/m3，那麼軟件自動計算物體的質量將會是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/m3. 那麼：
$$\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/m3
在工具→評估→質量屬性中能夠查看零件的慣性參數（若是單位或顯示的數值不合適能夠在選項中調整）
注意如下幾點：
- 重心/質心座標是相對於零件原點的
- SolidWorks列出了3個慣性張量，它們之間的區別就在於分別相對於不一樣的座標系：
- 相對於主軸座標系；其中的 Ix、Iy、Iz 三個向量表示主軸座標系相對於繪圖座標系（零件參考座標系）的姿態，即主軸座標系的 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]$$ - 相對於原點與主軸座標系重合，可是各軸與繪圖座標系一致的座標系
- 相對於繪圖座標系（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/m3，點擊按鈕進行自動計算：
經驗算，零件質量、質心位置以及慣性主矩與Solidworks的計算結果一致。區別是VREP中零件的固定參考座標系是默認在包圍盒的中心位置，目前還沒發現能夠更改這一參考點的地方。
參考：
Euler's Equations for Rotations in the Body-Fixed Frame
Introduction to Robotics - Mechanics and Control Chapter 6. Manipulator dynamics
Modern Robotics Mechanics, Planning, and Control Chapter 8. Dynamics of Open Chains
- 1. 牛頓迭代法(牛頓-拉弗森方法（Newton-Raphson method）)
- 2. 牛頓流體傳熱方程
- 3. MITE：歐拉歐拉歐拉歐拉(beta)
- 4. 歐拉四面體體積
- 5. 【梳理】離散數學 第15章 歐拉圖與哈密頓圖 15.1 歐拉圖 15.2 哈密頓圖
- 6. 牛頓方法
- 7. 牛頓法與擬牛頓法
- 8. 牛頓法解方程
- 9. 牛頓法與擬牛頓法學習筆記（一）牛頓法
- 10. 歐拉回路與歐拉路徑
- 更多相關文章...
- • 網站品質 教程 - 網站品質教程
- • ionic 下拉刷新 - ionic 教程
- • 常用的分佈式事務解決方案
- • Git五分鐘教程
-
每一个你不满意的现在，都有一个你没有努力的曾经。
- 1. 牛頓迭代法(牛頓-拉弗森方法（Newton-Raphson method）)
- 2. 牛頓流體傳熱方程
- 3. MITE：歐拉歐拉歐拉歐拉(beta)
- 4. 歐拉四面體體積
- 5. 【梳理】離散數學 第15章 歐拉圖與哈密頓圖 15.1 歐拉圖 15.2 哈密頓圖
- 6. 牛頓方法
- 7. 牛頓法與擬牛頓法
- 8. 牛頓法解方程
- 9. 牛頓法與擬牛頓法學習筆記（一）牛頓法
- 10. 歐拉回路與歐拉路徑