北京大学学报自然科学版   2016, Vol. 52 Issue(4): 708-712

文章信息

姚文莉
YAO Wenli
基于广义坐标形式的高斯最小拘束原理的多刚体系统动力学建模
Dynamical Modeling of Multi-Rigid-Body System Based on Gauss Principle of Least Constraint in Generalized Coordinates
北京大学学报(自然科学版), 2016, 52(4): 708-712
Acta Scientiarum Naturalium Universitatis Pekinensis, 2016, 52(4): 708-712

文章历史

收稿日期: 2015-11-23
修回日期: 2016-02-13
网络出版日期: 2016-07-12
基于广义坐标形式的高斯最小拘束原理的多刚体系统动力学建模
姚文莉     
青岛理工大学理学院, 青岛 266520
摘要: 通过采用动能及广义坐标显式的变分形式的高斯原理, 明确了广义坐标形式的高斯拘束中各项的含义, 以此建立以笛卡尔广义坐标表达的一般多刚体系统动力学问题的优化模型, 并研究利用上述模型列写其他坐标体系下的高斯拘束的方法。采用该方法可将多刚体系统的动力学问题变为求拘束极值的问题, 并且只要给出广义笛卡尔坐标与其他广义坐标之间的雅可比关系式, 便可方便地得到该坐标系统下的高斯拘束, 建模过程简单且具有更强的通用性。采用广义笛卡尔坐标及拉格朗日坐标, 对简单刚体的平面运动及定轴转动问题建立动力学优化模型, 并验证了该方法的有效性。
关键词: 高斯最小拘束原理     笛卡尔广义坐标     拉格朗日坐标体系     多刚体系统     动力学    
Dynamical Modeling of Multi-Rigid-Body System Based on Gauss Principle of Least Constraint in Generalized Coordinates
YAO Wenli     
School of Sciences, Qingdao Technological University, Qingdao 266520
Abstract: By variational Gauss principle in the form of kinetic energy and explicit express, the meaning of items in Gauss constraints is confirmed. The optimization model of dynamics of multi-rigid-body system in the form of Descartes generalized coordinates is established based on Gauss principle of least constraint in generalized coordinates. The method to write the Gauss constraints in other coordinate system according to the above model is studied. The dynamical problem of multi-rigid-body system can be turned into a problem for a minimum, and provided by the relation between the Descartes generalized coordinates and other coordinates system, the Gauss constraint in this kind of coordinates system is easy to be obtained. The modeling method is more convenient and is of universality. The dynamical models are established for planar motion and fixed axis rotation by Descartes generalized coordinate and Lagrange coordinate, and the validity of this method is proved.
Key words: Gauss principle of least constraint     Descartes generalized coordinates     Lagrange coordinate system     multi-rigid-body system     dynamics    

高斯最小拘束原理是多体系统动力学问题优化的重要原理之一, 它不需要建立系统的动力学方程, 而是以加速度为变量, 直接利用系统在每个时刻的坐标和速度得出运动规律。因其对约束形式没有要求, 所以可考察树系统、非树系统或单边约束的动力学问题, 还可将动力学分析与系统优化相结合, 解决带控制的多体系统动力学问题。将高斯最小拘束原理用于多刚体系统, 起源于机器人动力学研究[1]。目前, 利用高斯原理的多体系统建模方法均基于质点形式的高斯最小拘束原理, 需直接对多体系统内的每个物体单独列写拘束[2-5]。在用广义坐标表达的高斯最小拘束原理中, 对所选择的广义坐标没有要求, 所以建模过程简单且有更强的通用性, 但目前研究中, 广义坐标形式的高斯拘束原理中的符号含义不够清晰, 不适于程式化的建模[6-7]

本文采用动能及广义坐标显式的变分形式的高斯原理, 明确了广义坐标形式的高斯拘束中各项的含义, 并利用广义坐标形式的高斯最小拘束原理, 建立多刚体系统动力学的优化模型。在多刚体系统动力学问题的研究中, 针对不同的问题会采用不同的坐标系统, 与其他坐标系统相比, 广义笛卡尔坐标体系可以更方便地表达系统动力学量。因此, 本文研究利用上述模型列写其他坐标体系高斯拘束的方法, 根据上述方法, 只要能够得到其他坐标系统与广义笛卡尔坐标系统之间的关系表达式, 便可以方便地得到多刚体系统的拘束(该函数表达式中的变量均以矩阵形式给出)。

1 不同坐标形式的高斯最小拘束原理 1.1 以质点形式表达的高斯最小拘束原理

讨论由M个质点Pi(i=1, 2, …, M)组成的系统, 其拘束定义为$ Z = \frac{1}{2}\sum\limits_{i = 1}^M {{m_i}} {\left( {{{\ddot r}_i}-\frac{{{F_i}}}{{{m_i}}}} \right)^2} $, 将其展开为

$ \begin{array}{l} Z = \frac{1}{2}\sum\limits_{i = 1}^M {{m_i}} {\left( {{{\ddot r}_i}-\frac{{{F_i}}}{{{m_i}}}} \right)^2}\\ \;\;\; = \frac{1}{2}\sum\limits_{i = 1}^M {{m_i}} {{\ddot r}_i} \cdot {{\ddot r}_i}-\sum\limits_{i = 1}^M {{{\ddot r}_i} \cdot {F_i} + } \\ \;\;\;\;\;\;(与加速度的无关项), \end{array} $ (1)

则高斯最小拘束原理[8-9]表达为:在理想约束条件下, 系统在某瞬时, 在位置、速度和约束条件均相同但加速度不同的所有可能运动中, 真实运动使得拘束取极小值。

1.2 广义坐标形式的高斯最小拘束原理

研究受约束的系统, 其广义坐标可以表达为q1, q2, …, qn, 此组广义坐标可以是不独立的。将变分形式的高斯原理表达成Euler-Lagrange式的Gauss原理[10]

$ \sum\limits_{s = 1}^n {\left( {\frac{{\partial T}}{{\partial {q_s}}}-\frac{{\rm{d}}}{{{\rm{d}}t}}\frac{{\partial T}}{{\partial {{\dot q}_s}}} + {Q_s}} \right){\rm{\delta }}{{\ddot q}_s} = 0}, $ (2)

其中系统动能T

$ T = {T_2} + {T_1} + {T_0} = \frac{1}{2}\sum\limits_{s = 1}^n {\sum\limits_{k = 1}^n {{A_{ks}}{{\dot q}_s}} } {\dot q_k} + \sum\limits_{s = 1}^n {{B_s}} {\dot q_s} + {T_0}, $

则可将式(2)写为广义加速度的显式:

$ \sum\limits_{s = 1}^n {\left( {\sum\limits_{k = 1}^n {{A_{ks}}} {{\ddot q}_k} + {g_s}-{Q_s}} \right){\rm{\delta }}{{\ddot q}_s} = 0}。 $ (3)

定义下列矩阵:

$ \begin{array}{l} \mathit{\boldsymbol{A}} = {\left[{{A_{ks}}} \right]_{n \times n}}, \mathit{\boldsymbol{q}} = {[{q_{\rm{1}}}\;\;...\;\;{q_n}]^{\rm{T}}}, \;\\ \mathit{\boldsymbol{\dot q}} = {[{{\dot q}_{\rm{1}}}\;\;...\;\;{{\dot q}_n}]^{\rm{T}}}, \;\;\mathit{\boldsymbol{\ddot q}} = {[{{\ddot q}_{\rm{1}}}\;\;...\;\;{{\ddot q}_n}]^{\rm{T}}}, \end{array} $

$ \mathit{\boldsymbol{g}}(\mathit{\boldsymbol{q}},\;\mathit{\boldsymbol{\dot q}},\;t) = {[{g_{\rm{1}}}\;\;...\;\;{g_n}]^{\rm{T}}},\mathit{\boldsymbol{Q}} = {[{Q_{\rm{1}}}\;\;...\;\;{Q_n}]^{\rm{T}}}, $
$ {\mathit{\boldsymbol{Q}}_{\bf{s}}} = \sum\limits_{i = 1}^N {{F_i} \cdot \frac{{\partial {r_i}}}{{\partial {q_s}}}} $
$ \begin{array}{l} {g_s} = \sum\limits_{k = 1}^n {\sum\limits_{m = 1}^n {\left[{k, m;s} \right]{{\dot q}_k}} } {{\dot q}_m} -\sum\limits_{k = 1}^n {\left( {\frac{{\partial {B_k}}}{{\partial {q_s}}} -\frac{{\partial {B_s}}}{{\partial {q_k}}}} \right)\;} {{\dot q}_k} + \\ \;\;\;\;\;\;\frac{{\partial {B_s}}}{{\partial t}} -\frac{{\partial {T_0}}}{{\partial {q_s}}} + \sum\limits_{k = 1}^n {\frac{{\partial {A_{ks}}}}{{\partial t}}} {{\dot q}_k}, \end{array} $

其中,

$ [k, m;s] = \frac{1}{2}\left( {\frac{{\partial {A_{ks}}}}{{\partial {q_m}}} + \frac{{\partial {A_{ms}}}}{{\partial {q_k}}} -\frac{{\partial {A_{km}}}}{{\partial {q_s}}}} \right) $

为系数矩阵$ \mathit{\boldsymbol{A}} = {\left[{{A_{ks}}} \right]_{n \times n}} $的第一类Christoffel符号。

将式(3)写为矩阵形式:

$ {{\rm{(}}\mathit{\boldsymbol{A\ddot q}} + \mathit{\boldsymbol{g}}-\mathit{\boldsymbol{Q}})^{\rm{T}}}{\rm{ \mathsf{ δ} }}\mathit{\boldsymbol{\ddot q}} = 0。 $ (4)

假设在时刻t, q$ \mathit{\boldsymbol{\dot g}} $固定, 只有$ \mathit{\boldsymbol{\ddot g}} $变化, 故矩阵A, gQ均可看成不变量, 则由式(4)可以得到

$ G = \frac{1}{2}{{\rm{(}}\mathit{\boldsymbol{A\ddot q}} + \mathit{\boldsymbol{g}}-\mathit{\boldsymbol{Q}})^{\rm{T}}}{\mathit{\boldsymbol{A}}^{-1}}{\rm{(}}\mathit{\boldsymbol{A\ddot q}} + \mathit{\boldsymbol{g}}-\mathit{\boldsymbol{Q}})。 $ (5)

因此, 高斯最小拘束可以取为

$ {\rm{ \mathsf{ δ} }}[{{\rm{(}}\mathit{\boldsymbol{A\ddot q}} + \mathit{\boldsymbol{g}}-\mathit{\boldsymbol{Q}})^{\rm{T}}}{\mathit{\boldsymbol{A}}^{-1}}{\rm{(}}\mathit{\boldsymbol{A\ddot q}} + \mathit{\boldsymbol{g}}-\mathit{\boldsymbol{Q}})] = 0。 $

即在时刻t, q$ \mathit{\boldsymbol{\dot g}} $固定的情况下, 在所有满足约束的可能运动中, 实际运动所对应的高斯最小拘束最小。

在上述定理的推导中, 广义加速度并未要求是独立的, 适用性很广, 对完整系统、非完整系统以及单、双边约束都成立。

2 采用笛卡尔广义坐标形式的多刚体系统高斯最小拘束原理

设多刚体系统由N个刚体Bi组成, 地球为零刚体B0。取定一个惯性参考基e(0)和每个刚体的连体基e(i), e(i)的原点与质心Ci重合。为了确定系统内每个刚体相对于惯性基的位形, 可以用质心Ci的位置矢径的3个分量(xc, yc, zc)i确定位置, 用连体基的3个欧拉角(ψ, θ, φ)i确定方位。

将3个平动坐标和3个转动坐标写成6×1矢量列阵:

$ \mathit{\boldsymbol{q}} = {[\mathit{\boldsymbol{q}}_{\rm{1}}^{\rm{T}}\;\;...\;\;q_N^{\rm{T}}]^{\rm{T}}}, $

根据欧拉运动学方程:

$ {\left( {\begin{array}{*{20}{c}} {{\omega _{x'}}}\\ {{\omega _{y'}}}\\ {{\omega _{z'}}} \end{array}} \right)_i} = {\mathit{\boldsymbol{B}}_i}({\theta _i}, {\phi _i})\left( {\begin{array}{*{20}{c}} {{{\dot \psi }_i}}\\ {{{\dot \theta }_i}}\\ {{{\dot \phi }_i}} \end{array}} \right), $

其中,

$ {\mathit{\boldsymbol{B}}_i}({\theta _i}, \;{\phi _i}) = {\left[{\begin{array}{*{20}{c}} {\sin \theta \sin \phi }&{\cos \phi }&0\\ {\sin \theta \cos \phi }&{-\sin \phi }&0\\ {\cos \theta }&0&1 \end{array}} \right]_i}, $

系统动能可写为

$ \begin{array}{l} {T_i} = \frac{1}{2}{[\dot x, \dot y, \dot z]_i}{\left[{\begin{array}{*{20}{c}} m&0&0\\ 0&m&0\\ 0&0&m \end{array}} \right]_i}{\left[{\begin{array}{*{20}{c}} {\dot x}\\ {\dot y}\\ {\dot z} \end{array}} \right]_i} + \\ \;\;\;\;\;\frac{1}{z}{\left[{\begin{array}{*{20}{c}} {{\omega _{x'}}}&{{\omega _{y'}}}&{{\omega _{z'}}} \end{array}} \right]_i}{\left[{\begin{array}{*{20}{c}} {{J_{x'}}}&0&0\\ 0&{{J_{y'}}}&0\\ 0&0&{{J_{z'}}} \end{array}} \right]_i}{\left[{\begin{array}{*{20}{c}} {{\omega _{x'}}}\\ {{\omega _{y'}}}\\ {{\omega _{z'}}} \end{array}} \right]_i}\\ \;\; = \frac{1}{2}\mathit{\boldsymbol{\dot q}}_i^{\rm{T}}{\mathit{\boldsymbol{M}}_i}{{\mathit{\boldsymbol{\dot q}}}_i}\;。 \end{array} $

$ {{\mathit{\boldsymbol{J'}}}_i} = {\left[{\begin{array}{*{20}{c}} {{J_{x'}}}&0&0\\ 0&{{J_{y'}}}&0\\ 0&0&{{J_{z'}}} \end{array}} \right]_i}{\rm{, }}\;{\mathit{\boldsymbol{m}}_i} = {\left[{\begin{array}{*{20}{c}} m&0&0\\ 0&m&0\\ 0&0&m \end{array}} \right]_i}, $, 定义惯性矩阵$ {\mathit{\boldsymbol{M}}_i} = \left[{\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{m}}_i}}&0\\ 0&{{\mathit{\boldsymbol{J}}_i}} \end{array}} \right] $, 其中$ {\mathit{\boldsymbol{J}}_i} = {\mathit{\boldsymbol{B}}_i}^{\rm{T}}\mathit{\boldsymbol{J}}{\mathit{\boldsymbol{'}}_i}{\mathit{\boldsymbol{B}}_i} $, 可以证明, 只要J'对称, 则$ {\mathit{\boldsymbol{J}}_i} = {\mathit{\boldsymbol{B}}_i}^{\rm{T}}\mathit{\boldsymbol{J}}{\mathit{\boldsymbol{'}}_i}{\mathit{\boldsymbol{B}}_i} $对称。因为在广义笛卡尔坐标系下, J'总是对称的, 故$ {\mathit{\boldsymbol{J}}_i} = {\mathit{\boldsymbol{B}}_i}^{\rm{T}}\mathit{\boldsymbol{J}}{\mathit{\boldsymbol{'}}_i}{\mathit{\boldsymbol{B}}_i} $为对称矩阵。

N个刚体组成的多刚体位形由6N个广义坐标确定, 可以写成6N×1的位置矢量列阵:

$ \mathit{\boldsymbol{q}} = {[\mathit{\boldsymbol{q}}_{\rm{1}}^{\rm{T}}\;\;...\;\;q_N^{\rm{T}}]^{\rm{T}}}, $
$ \mathit{\boldsymbol{\dot X}} = \left[{\begin{array}{*{20}{c}} {{{\mathit{\boldsymbol{\dot x}}}_{\bf{1}}}}\\ {{{\mathit{\boldsymbol{\dot x}}}_{\bf{2}}}}\\ \vdots \\ {{{\mathit{\boldsymbol{\dot x}}}_N}} \end{array}} \right], {\rm{ }}{\bf{ }}\mathit{\boldsymbol{M}} = {\rm{diag}}({\mathit{\boldsymbol{M}}_1}, \;{\mathit{\boldsymbol{M}}_2} \cdots, \;{\mathit{\boldsymbol{M}}_N}), $
$ \mathit{\boldsymbol{Q}} = {[\mathit{\boldsymbol{Q}}_{\rm{1}}^{\rm{T}}\;\;...\;\;\mathit{\boldsymbol{Q}}_N^{\rm{T}}]^{\rm{T}}}, $

其中相对于刚体i的广义笛卡尔坐标的广义力列阵[11]

$ {\mathit{\boldsymbol{Q}}_i} = \left[{{F_x}, {F_y}, {F_y}, {M_z}, {M_N}, {M_{z'}}} \right]_i^{\rm{T}}, $ (6)

列阵内各元素分别为外力主矢在定坐标系中的投影及外力对3根转动轴的主矩。

针对广义笛卡尔坐标形式的高斯最小拘束为

$ \begin{array}{l} G = \frac{1}{2}{{\rm{(}}\mathit{\boldsymbol{M\ddot X}} + \mathit{\boldsymbol{g}} - \mathit{\boldsymbol{Q}})^{\rm{T}}}{\mathit{\boldsymbol{M}}^{ - 1}}{\rm{(}}\mathit{\boldsymbol{M\ddot X}} + \mathit{\boldsymbol{g}} - \mathit{\boldsymbol{Q}})\\ \;\;\; = \frac{1}{2}{{\ddot X}^{\rm{T}}}\mathit{\boldsymbol{M}}\ddot X + {(\mathit{\boldsymbol{g}} - \mathit{\boldsymbol{Q}})^{\rm{T}}}\ddot X + (与加速度无关项), \end{array} $

故可取高斯最小拘束函数为

$ G = \frac{1}{2}{{\mathit{\boldsymbol{\ddot X}}}^{\rm{T}}}\mathit{\boldsymbol{M\ddot X}} + {(\mathit{\boldsymbol{g}}-\mathit{\boldsymbol{Q}})^{\rm{T}}}\mathit{\boldsymbol{\ddot X}}。 $ (7)

因此, 高斯最小拘束原理可表达为:在时刻t, q$ \mathit{\boldsymbol{\dot g}} $固定的情况下, 在所有满足约束(对约束类型没有限制, 可以为完整或非完整及单双边约束)的可能运动中, 实际运动高斯拘束取得最小值。

以广义笛卡尔坐标表达的优点在于其惯性矩阵规范简单, 但在实际应用中, 也常常需要建立其他坐标系统的动力学模型。例如, 对于自动机械、机器人等带控制的系统, 反馈控制变量是系统内刚体间连接铰链的相对运动变量, 但在用广义笛卡尔坐标建立的动力学模型中不出现该变量, 因此这类模型难以应用于带控制的系统。下面, 讨论如何根据现有广义笛卡尔坐标下的高斯拘束得到其他广义坐标系统下高斯拘束的方法。

3 其他广义坐标系统下的高斯最小拘束

N个刚体组成的多刚体系统受到连接铰链产生的约束, 这6N个广义笛卡尔坐标一般不是独立的, 可以选另一套广义坐标$ \mathit{\boldsymbol{q'}} = {[{{q'}_1}, \;\;{{q'}_2}, \;..., \;{{q'}_n}]^{\rm{T}}} $作为广义坐标列阵(如选择系统内刚体之间连接铰链的相对运动变量, 即铰链坐标)。如果多刚体系统是开链系统, 则铰链坐标是独立的, 若是闭链系统, 则坐标不独立, 坐标数目等于它的开链缩减系统的自由度数目, 约束条件由回路闭合条件决定。

若在广义笛卡尔坐标和该广义坐标之间存在下列关系:

$ \mathit{\boldsymbol{X}} = \mathit{\boldsymbol{X}}(\mathit{\boldsymbol{q'}}, \;t), \;{x_i} = {x_i}({{q'}_1}, \;..., \;{{q'}_n}, \;t), $

则下式成立:

$ \left\{ \begin{array}{l} \mathit{\boldsymbol{\dot X}} = \frac{{\partial \mathit{\boldsymbol{X}}}}{{\partial \mathit{\boldsymbol{q'}}}}\mathit{\boldsymbol{\dot q'}} + \frac{{\partial \mathit{\boldsymbol{X}}}}{{\partial t}}\\ \mathit{\boldsymbol{\ddot X}} = \frac{{\partial \mathit{\boldsymbol{X}}}}{{\partial \mathit{\boldsymbol{q'}}}}\mathit{\boldsymbol{\ddot q'}} + \frac{{{\partial ^2}\mathit{\boldsymbol{X}}}}{{\partial {{\mathit{\boldsymbol{q'}}}^2}}}\mathit{\boldsymbol{\dot q'}} + \frac{{{\partial ^2}\mathit{\boldsymbol{X}}}}{{\partial \mathit{\boldsymbol{\dot q}}\partial t}}\mathit{\boldsymbol{\dot q'}} + \frac{{{\partial ^2}\mathit{\boldsymbol{X}}}}{{\partial \mathit{\boldsymbol{q'}}\partial t}} + \frac{{{\partial ^2}\mathit{\boldsymbol{X}}}}{{\partial {t^2}}}, \\ {{\dot x}_i} = \sum\limits_{l = 1}^n {\frac{{\partial {x_i}}}{{\partial {{q'}_l}}}{{\dot q}_l} + \frac{{\partial {x_i}}}{{\partial t}}} \\ {{\ddot x}_i} = \sum\limits_{l = 1}^n {\frac{{\partial {x_i}}}{{\partial {{q'}_l}}}{{\ddot q'}_l}} + \sum\limits_{l = 1}^n {\frac{{{\partial ^2}{x_i}}}{{\partial q_l^{'2}}}{{\dot q'}_l}} + \sum\limits_{l = 1}^n {\frac{{{\partial ^2}x}}{{\partial {{q'}_l}\partial t}}{{\dot q'}_l}} + \\ \;\;\;\;\;\;\sum\limits_{l = 1}^n {\frac{{{\partial ^2}x}}{{\partial {{q'}_l}\partial t}}} + \frac{{{\partial ^2}x}}{{\partial {t^2}}}, \end{array} \right. $ (8)
$ {\rm{ \mathsf{ δ} }}\mathit{\boldsymbol{\ddot X}} = \frac{{\partial \mathit{\boldsymbol{X}}}}{{\partial \mathit{\boldsymbol{q'}}}}{\rm{ \mathsf{ δ} }}\mathit{\boldsymbol{\ddot q'}}。 $ (9)

根据式(4), 可得到下列变分表达式:

$ {{\rm{(}}\mathit{\boldsymbol{M\ddot X}} + \mathit{\boldsymbol{g}}-\mathit{\boldsymbol{Q}})^{\rm{T}}}{\rm{ \mathsf{ δ} }}\mathit{\boldsymbol{\ddot X}} = {\rm{0}}。 $

将式(8)和(9)代入上式得

$ {\left[{\mathit{\boldsymbol{M}}\frac{{\partial \mathit{\boldsymbol{X}}}}{{\partial \mathit{\boldsymbol{q'}}}}\ddot q' + \mathit{\boldsymbol{H}}-\mathit{\boldsymbol{Q}}} \right]^{\rm{T}}}\frac{{\partial \mathit{\boldsymbol{X}}}}{{\partial \mathit{\boldsymbol{q'}}}}{\rm{ \mathsf{ δ} }}\mathit{\boldsymbol{\ddot q'}} = 0, $ (10)

其中,

$ \mathit{\boldsymbol{H}} = \mathit{\boldsymbol{M}}\frac{{{\partial ^2}\mathit{\boldsymbol{X}}}}{{\partial {{\mathit{\boldsymbol{q'}}}^2}}}\mathit{\boldsymbol{\dot q'}} + \mathit{\boldsymbol{M}}\frac{{{\partial ^2}\mathit{\boldsymbol{X}}}}{{\partial \mathit{\boldsymbol{q'}}\partial t}}\mathit{\boldsymbol{\dot q'}} + \mathit{\boldsymbol{M}}\frac{{{\partial ^2}\mathit{\boldsymbol{X}}}}{{\partial \mathit{\boldsymbol{q'}}\partial t}} + \mathit{\boldsymbol{M}}\frac{{{\partial ^2}\mathit{\boldsymbol{X}}}}{{\partial {t^2}}} + g, $

式(10)可化为

$ \begin{array}{l} \left[{{{\mathit{\boldsymbol{\ddot q'}}}^{\rm{T}}}{{\left( {\frac{{\partial \mathit{\boldsymbol{X}}}}{{\partial \mathit{\boldsymbol{q'}}}}} \right)}^{\rm{T}}}{\mathit{\boldsymbol{M}}^{\rm{T}}}\frac{{\partial \mathit{\boldsymbol{X}}}}{{\partial \mathit{\boldsymbol{q'}}}} + {\mathit{\boldsymbol{H}}^{\rm{T}}}\frac{{\partial \mathit{\boldsymbol{X}}}}{{\partial \mathit{\boldsymbol{q'}}}}-{\mathit{\boldsymbol{Q}}^{\rm{T}}}\frac{{\partial \mathit{\boldsymbol{X}}}}{{\partial \mathit{\boldsymbol{q'}}}}} \right]\delta \mathit{\boldsymbol{\ddot q'}} = 0, \\ {\left[{{{\left( {\frac{{\partial \mathit{\boldsymbol{X}}}}{{\partial \mathit{\boldsymbol{q'}}}}} \right)}^{\rm{T}}}{\mathit{\boldsymbol{M}}^{\rm{T}}}\frac{{\partial \mathit{\boldsymbol{X}}}}{{\partial \mathit{\boldsymbol{q'}}}}\mathit{\boldsymbol{\ddot q'}} + {{\left( {\frac{{\partial \mathit{\boldsymbol{X}}}}{{\partial \mathit{\boldsymbol{q'}}}}} \right)}^{\rm{T}}}\mathit{\boldsymbol{H}}-{{\left( {\frac{{\partial \mathit{\boldsymbol{X}}}}{{\partial \mathit{\boldsymbol{q'}}}}} \right)}^{\rm{T}}}\mathit{\boldsymbol{Q}}} \right]^{\rm{T}}}\delta \mathit{\boldsymbol{\ddot q'}} = 0。\; \end{array} $

$ \mathit{\boldsymbol{M'}} = {\left( {\frac{{\partial \mathit{\boldsymbol{X}}}}{{\partial \mathit{\boldsymbol{q'}}}}} \right)^{\rm{T}}}{\mathit{\boldsymbol{M}}^{\rm{T}}}\frac{{\partial \mathit{\boldsymbol{X}}}}{{\partial \mathit{\boldsymbol{q'}}}} = {{\mathit{\boldsymbol{M'}}}^{\rm{T}}} $(对称矩阵), $ \mathit{\boldsymbol{g'}} = {\left( {\frac{{\partial \mathit{\boldsymbol{X}}}}{{\partial \mathit{\boldsymbol{q'}}}}} \right)^{\rm{T}}}\mathit{\boldsymbol{H}}, \mathit{\boldsymbol{Q'}} = {\left( {\frac{{\partial \mathit{\boldsymbol{X}}}}{{\partial \mathit{\boldsymbol{q}}}}} \right)^{\rm{T}}}\mathit{\boldsymbol{Q}} $, 则该广义坐标形式的高斯最小拘束:

$ \begin{array}{l} G = \frac{1}{2}{\left[{{{\left( {\frac{{\partial \mathit{\boldsymbol{X}}}}{{\partial \mathit{\boldsymbol{q'}}}}} \right)}^{\rm{T}}}{\mathit{\boldsymbol{M}}^{\rm{T}}}\frac{{\partial \mathit{\boldsymbol{X}}}}{{\partial \mathit{\boldsymbol{q'}}}}\mathit{\boldsymbol{\ddot q'}} + {{\left( {\frac{{\partial \mathit{\boldsymbol{X}}}}{{\partial \mathit{\boldsymbol{q'}}}}} \right)}^{\rm{T}}}\mathit{\boldsymbol{H}}-{{\left( {\frac{{\partial \mathit{\boldsymbol{X}}}}{{\partial \mathit{\boldsymbol{q'}}}}} \right)}^{\bf{T}}}\mathit{\boldsymbol{Q}}} \right]^{\rm{T}}} \cdot \\ {\left[{{{\left( {\frac{{\partial \mathit{\boldsymbol{X}}}}{{\partial \mathit{\boldsymbol{q'}}}}} \right)}^{\bf{T}}}{\mathit{\boldsymbol{M}}^{\rm{T}}}\frac{{\partial \mathit{\boldsymbol{X}}}}{{\partial \mathit{\boldsymbol{q'}}}}} \right]^{ - 1}}\left[{{{\left( {\frac{{\partial \mathit{\boldsymbol{X}}}}{{\partial \mathit{\boldsymbol{q'}}}}} \right)}^{\rm{T}}}{\mathit{\boldsymbol{M}}^{\rm{T}}}\frac{{\partial \mathit{\boldsymbol{X}}}}{{\partial \mathit{\boldsymbol{q'}}}}\mathit{\boldsymbol{\ddot q'}} + } \right.\\ \left. {\;\;\;{{\left( {\frac{{\partial \mathit{\boldsymbol{X}}}}{{\partial \mathit{\boldsymbol{q'}}}}} \right)}^{\rm{T}}}\mathit{\boldsymbol{H}}-{{\left( {\frac{{\partial \mathit{\boldsymbol{X}}}}{{\partial \mathit{\boldsymbol{q'}}}}} \right)}^{\rm{T}}}} \right]{\rm{ }}\\ {\rm{ = }}\frac{1}{2}{{\rm{(}}\mathit{\boldsymbol{M'\ddot q'}} + \mathit{\boldsymbol{g'}} -\mathit{\boldsymbol{Q'}})^{\rm{T}}}{{\mathit{\boldsymbol{M'}}}^{ -1}}{\rm{(}}\mathit{\boldsymbol{M'\ddot q'}} + \mathit{\boldsymbol{g'}} -\mathit{\boldsymbol{Q'}})\\ = \frac{1}{2}{{\mathit{\boldsymbol{\ddot q'}}}^{\rm{T}}}{{\mathit{\boldsymbol{M'}}}^{\rm{T}}}\mathit{\boldsymbol{\ddot q'}} + {{\rm{(}}\mathit{\boldsymbol{g'}} - \mathit{\boldsymbol{Q'}})^{\rm{T}}}\mathit{\boldsymbol{\ddot q'}} + \left( {与加速度无关项}\right), \end{array} $

故取高斯拘束为

$ G = \frac{1}{2}{{\mathit{\boldsymbol{\ddot q'}}}^{\rm{T}}}{{\mathit{\boldsymbol{M'}}}^{\rm{T}}}\mathit{\boldsymbol{\ddot q'}} + {{\rm{(}}\mathit{\boldsymbol{g'}}-\mathit{\boldsymbol{Q'}})^{\rm{T}}}\mathit{\boldsymbol{\ddot q'}}。 $ (11)

因此, 只要给出广义笛卡尔坐标与其他坐标之间的雅可比关系式, 便可利用式(11), 方便地得到用其他广义坐标表达的高斯最小拘束。

4 例子 4.1 计算平面运动刚体的高斯拘束

取广义笛卡尔坐标系统$ \mathit{\boldsymbol{X}} = {[{x_c}, \;\;{y_c}, \;\;{z_c}, \;\;\psi, \;\;\theta, \;\varphi]^{\rm{T}}}, $对于平面刚体的问题, 总有$ {z_c} \equiv 0, \psi \equiv 0, \theta \equiv 0 $成立, 则

$ \begin{array}{l} {\mathit{\boldsymbol{B}}_i}({\theta _i}, \;{\phi _i}) = \left[{\begin{array}{*{20}{c}} {\sin \theta \sin \phi }&{\cos \phi }&0\\ {\sin \theta \cos \phi }&{-\sin \phi }&0\\ {\cos \theta }&0&1 \end{array}} \right] = \left[{\begin{array}{*{20}{c}} 0&{\cos \phi }&0\\ 0&{-\sin \phi }&0\\ 1&0&1 \end{array}} \right], \\ \mathit{\boldsymbol{g}} = {\bf{0}}, \end{array} $

其高斯最小拘束为

$ \begin{array}{l} G = \frac{1}{2}{{\mathit{\boldsymbol{\ddot X}}}^{\rm{T}}}\mathit{\boldsymbol{M\ddot X}} + {(\mathit{\boldsymbol{g}}-\mathit{\boldsymbol{Q}})^{\rm{T}}}\mathit{\boldsymbol{\ddot X}}\\ {\rm{ }} = \frac{1}{2}m(\ddot x_c^2 + \ddot y_c^2) + \frac{1}{2}{J_z}{{\dot \phi }^2}-{F_x}{{\ddot x}_c}-{F_y}{{\ddot y}_c} - {M_z}\ddot \phi 。 \end{array} $

高斯最小拘束满足驻值条件$ {\rm{\delta }}G = 0 $, 即$ \frac{{\partial G}}{{\partial \ddot x}} = 0, $$ \frac{{\partial G}}{{\partial \ddot y}} = 0, \; $$ \frac{{\partial G}}{{\partial \ddot \varphi }} = 0 $, 则下式成立:

$ \left\{ \begin{array}{l} m{{\ddot x}_c} = {F_x}, \\ m{{\ddot y}_c} = {F_y}, \\ {J_z}\ddot \phi = {M_z}。\; \end{array} \right. $

上式与刚体平面微分方程[12]一致。

4.2 计算定轴转动刚体的高斯拘束

平面运动刚体在运动过程中存在定点A, 则系统取广义坐标q'=φ, 该广义坐标与广义笛卡尔坐标之间的关系为

$ {x_c} = l \cdot \sin q', \;{y_c} = l \cdot \cos q', \;\;\varphi = q'。 $

计算高斯拘束(11)中的各项:

$ \mathit{\boldsymbol{M'}} = {\left( {\frac{{\partial \mathit{\boldsymbol{X}}}}{{\partial \mathit{\boldsymbol{q'}}}}} \right)^{\rm{T}}}{\mathit{\boldsymbol{M}}^{\rm{T}}}\frac{{\partial \mathit{\boldsymbol{X}}}}{{\partial \mathit{\boldsymbol{q'}}}} = [m \cdot {l^2} + {J_C}] = [{J_A}], $
$ \mathit{\boldsymbol{g'}} = {\left( {\frac{{\partial \mathit{\boldsymbol{X}}}}{{\partial \mathit{\boldsymbol{g'}}}}} \right)^{\rm{T}}}\mathit{\boldsymbol{H}} = {\bf{0}}, $
$ {{\mathit{\boldsymbol{Q'}}}^{\rm{T}}} = {\mathit{\boldsymbol{Q}}^{\rm{T}}}\left( {\frac{{\partial \mathit{\boldsymbol{X}}}}{{\partial \mathit{\boldsymbol{q'}}}}} \right) = [{F_x}l \cdot \cos \theta-{F_y}l \cdot \sin \theta + {M_C}] = [{M_A}], $

则在此广义坐标下的高斯拘束为

$ G = \frac{1}{2}{{\mathit{\boldsymbol{\ddot q'}}}^{\rm{T}}}{{\mathit{\boldsymbol{M'}}}^{\rm{T}}}\mathit{\boldsymbol{\ddot q'}} + {(\mathit{\boldsymbol{g'}} - \mathit{\boldsymbol{Q'}})^{\rm{T}}}\mathit{\boldsymbol{\ddot q'}} = \frac{1}{2}{J_A}{{\ddot \theta }^2} - {M_A}\ddot \theta 。 $

根据δG=0, 即$ \frac{{\partial G}}{{\partial \mathit{\boldsymbol{\ddot q'}}}} = 0 $, 可以得到, $ {J_A}\ddot \theta = {M_A} $为刚体定轴转动的运动微分方程[11]

5 结语

高斯最小拘束原理是研究多体系统动力学问题优化方法的重要原理之一。本文采用动能及广义加速度显式的变分形式的高斯原理, 明确了广义坐标形式的高斯最小拘束原理高斯拘束中各项的含义, 并推导了广义笛卡尔坐标表达的一般多刚体系统动力学问题的优化模型, 建立了各种广义坐标与广义笛卡尔坐标表达的高斯拘束之间的关系。当根据需要选取不同的广义坐标时, 可以通过优化方法解决多刚体系统动力学问题建模, 过程简单且规范。在刚体平面运动及刚体做定轴转动的例子中, 计算了两种坐标系下的高斯拘束, 并通过高斯最小拘束原理的必要条件验证了本文的优化模型, 它与典型的动力学微分方程一致, 间接证明了该方法的有效性。

参考文献
[1] 袁士杰, 吕哲勤. 多刚体系统动力学. 北京: 北京理工大学出版社, 1991 .
[2] 刘延柱. 基于高斯原理的多体系统动力学建模. 力学学报 , 2014, 46 (6) : 940–945.
[3] 董龙雷, 闫桂荣, 杜彦亭, 等. 高斯最小拘束原理在一类刚柔耦合系统分析中的应用. 兵工学报 , 2001, 22 (3) : 347–351.
[4] 史跃东, 王德石. 舰炮振动的高斯最小拘束分析方法. 海军工程大学学报 , 2009, 21 (5) : 1–5.
[5] 郝名望, 叶正寅. 单柔体与多柔体动力学的高斯最小拘束原理. 广西大学学报 , 2011, 36 (2) : 195–204.
[6] 姚文莉, 戴葆青. 广义坐标形式的高斯最小拘束原理及其推广. 力学与实践 , 2014, 36 (6) : 779–785.
[7] Kalaba R E, Udwadia F E. Equations of motion for nonholonomic, constrained dynamical systems via Gauss's Principle. J Appl Mech , 1993, 60 (3) : 662–668 DOI:10.1115/1.2900855 .
[8] 陈滨. 分析动力学. 2版. 北京: 北京大学出版社, 2012 .
[9] 刘延柱. 高等动力学. 北京: 高等教育出版社, 2000 .
[10] 梅凤翔, 刘端, 罗勇. 高等分析力学. 北京: 北京理工大学出版社, 1991 .
[11] 黄昭度, 钟奉俄. 工程系统分析力学. 北京: 高等教育出版社, 1992 .
[12] 哈尔滨工业大学理论力学教研室. 理论力学. 7版. 北京: 高等教育出版社, 2009 .