CN105519267B - 一种变构型航天器多刚体动力学建模方法 - Google Patents

Abstract

本发明公开了一种变构型航天器多刚体动力学建模方法，根据变构型航天器的拓扑构型分析，建立变构型航天器的拓扑结构模型，建立坐标系，对变构型航天器进行运动学分析，根据单个变构型航天器舱段的运动学参量，建立邻接舱段运动学递推关系，描述变构型航天器运动，根据已知运动规律的舱段和步骤二中得到的邻接舱段递推关系，建立变构型航天器系统动力学方程，对变构型航天器的系统动力学方程积分，以变构型航天器所受外力作为已知约束条件，求得描述变构型航天器的动力学参量，建立变构型航天器的多刚体动力学模型。本发明能减少人工拓扑分析，便于计算机编程实现拓扑构型分析。<pb pnum="1" />

Description

一种变构型航天器多刚体动力学建模方法

技术领域

本发明属于航天器动力学领域，涉及一种变构型航天器的拓扑构型描述及动力学建模方法。

背景技术

随着航天任务的复杂化、多样化及突发性，对航天器提出了低成本、多任务、高机动的要求，变构型航天器成为未来航天器的主流发展方向。

变构型航天器是指在利用大规模舱段模块堆栈设计，并能够在空间展开并自重构的多舱段航天器。它是一种高维多自由度动力学系统，由于热环境、舱段模块的同时转动、以及控制力的加入，形成多场强耦合，加剧了变构型航天器姿态控制的难度。在处理这类系统的建模问题中，总的来说有两类方法：一是柔性多体动力学建模方法，利用有限元划分网格，模态叠加及模态截断方法，在建模过程中引入柔性耦合项，该方法适合于大变形下的多体动力学问题；另一种，则是多刚体动力学建模方法，采用图论方法描述变构型航天器拓扑结构，将柔性振动作为扰动项与刚体模型叠加。前一种方法虽然考虑了刚柔耦合效应，但是有限元离散和模态叠加使得计算量庞大，对星载计算机的要求较高，实时性较差。后者采用刚体假设，针对一般规则结构的空间柔性航天器来说结果足够准确。再者，多刚体动力学已经历十余年发展，具有很成熟的数学工具支持，便于编程处理，实现自动化分析，相比之下，能更好的满足变构型航天器实时控制的要求。

针对后一种方法，变构型航天器的拓扑构型描述目前仍然需要地面人工完成，每次构型发生变化后，都需要重新进行工作量巨大的拓扑分析。这将严重的限制航天器的自主性、快速性、实时控制等在轨操作性能。

总的来说，现有的变构型航天器柔性多体动力学建模方法中，对减少人工拓扑分析工作问题处理还没有明确的方法和解决方案。

发明内容

为了克服现有技术的不足，本发明提供一种变构型航天器多刚体动力学建模方法，能减少人工拓扑分析，便于计算机编程实现拓扑构型分析。

本发明解决其技术问题所采用的技术方案包括以下步骤：

步骤一、根据变构型航天器的拓扑构型分析，建立变构型航天器的拓扑结构模型，获得变构型航天器的邻接矩阵D、关联矩阵S、内接物体矩阵L、内接物体数组L；

步骤二、建立坐标系，对变构型航天器进行运动学分析，包括以下步骤：

步骤A、建立参考惯性坐标系、体坐标系、铰坐标系，

建立右手惯性坐标系OXYZ，记作参考基e，基点O取在航天器初始构型D₀的整体质心处，X轴、Y轴、Z轴分别沿初始构型下的柯西惯性椭球长、中、短轴方向，满足右手坐标系定义；邻接两舱段中，舱段B_j为舱段B_i的内接舱段，即j＝L(i)：两舱段由铰H_i相连，在动力学模型中，忽略铰的形状和尺寸，简化为一个几何点，称为铰点，P与Q分别为铰点在舱段B_i与B_j上的位置；以舱段B_i与B_i的质心C_i与C_j为基点分别建立连体基，即舱段B_i与B_j的体坐标系，分别记作e ⁱ和e ^j，它们关于e的方向余弦阵分别记为A ⁱ与A ^j；建立以Q为原点的铰H_i的铰坐标系，即铰H_i的本地基，固结于点Q，记为过P建立铰H_i的动基e ^h，它固结于点P；

步骤B、采用广义坐标描述单舱段运动学参量，

定义广义坐标q _i(i＝1，…，Num)，组成列阵Num为铰的个数，根据规则标号法，对于树状拓扑结构的航天器系统，铰个数Num等于舱段个数N，刚体舱段B_i相对内接舱段B_j的相对位移相对速度转动角速度其中H ^′hT、H ^′ΩT是单位矢量动基e ^h的本地坐标阵，仅和铰的形式有关，统一用上标“′”表示相对于当地坐标系(即B_i浮动体坐标系)的矢量坐标阵，无上标表示相对于绝对坐标系的矢量坐标阵；

步骤三、根据单个变构型航天器舱段的运动学参量，建立邻接舱段运动学递推关系，描述变构型航天器运动，包括以下内容；

变构型航天器系统内各舱段角速度组成的列阵ω＝(ω₁，…，ω_N)^T，则有

$\underset{\&OverBar;}{\mathrm{\&omega;}}=\underset{\&OverBar;}{\mathrm{\&beta;}}\underset{\&OverBar;}{\stackrel{\&CenterDot;}{q}}+{\mathrm{\&omega;}}_0{\underset{\&OverBar;}{l}}_N,\underset{\&OverBar;}{\stackrel{\&CenterDot;}{\mathrm{\&omega;}}}=\underset{\&OverBar;}{\mathrm{\&beta;}}\underset{\&OverBar;}{\stackrel{\&CenterDot;\&CenterDot;}{q}}+\underset{\&OverBar;}{\mathrm{\&sigma;}}$

其中

各舱段质心位置列阵r＝(r₁，…，r_N)^T＝d ^*T 1 _N+(r₀+h₁)1 _N，质心速度其中h₁为六自由度虚铰H₁的滑移铰矢量，α＝-(H ^Ω T×d ^*+H ^h T)^T，

$\underset{\&OverBar;}{\mathrm{\&upsi;}}=-{\underset{\&OverBar;}{d}}^{*T}\&times;{\mathrm{\&omega;}}_0{\underset{\&OverBar;}{l}}_N+({\mathrm{\&omega;}}_0\&times;h_1+{\stackrel{\&CenterDot;}{r}}_0){\underset{\&OverBar;}{l}}_N,$

d^*为广义通路矢量，其中c_ii为舱段B_i质心C_i指向B_i上铰点P的矢量；

步骤四、根据已知运动规律的舱段和步骤二中得到的邻接舱段递推关系，建立变构型航天器系统动力学方程，包括以下内容：

根据速度变分形式的动力学普遍方程，可得微分形式的系统动力学方程：

$\underset{\&OverBar;}{Z}\underset{\&OverBar;}{\stackrel{\&CenterDot;\&CenterDot;}{q}}=\underset{\&OverBar;}{z}$

其中Z为广义质量矩阵，z为广义力阵：

Z＝α ^*T mα ^*+β ^T Jβ

z＝α ^*T·(μF ^Out-mμw)+β ^T·(M ^O-J·σ-ε)+F ^eq+F ^τq

式中F ^Out为外力约束己知条件，F ^τq为直接作用于铰点的约束力元，F ^eq为作用于舱段上铰点以外的非约束力元，J为中心惯量张量J _i组成的对角阵，

$\underset{\&OverBar;}{\mathrm{\&mu;}}={\underset{\&OverBar;}{l}}_N-\frac{1}{m_s}\underset{\&OverBar;}{m}{\underset{\&OverBar;}{l}}_N{\underset{\&OverBar;}{l}}_N^T,{\underset{\&OverBar;}{\mathrm{\&alpha;}}}^{*}={\underset{\&OverBar;}{\mathrm{\&mu;}}}^T\underset{\&OverBar;}{\mathrm{\&alpha;}},{\underset{\&OverBar;}{\mathrm{\&mu;}}}^T\underset{\&OverBar;}{w}=-{\left({\underset{\&OverBar;}{d}}^{*}\underset{\&OverBar;}{\mathrm{\&mu;}}\right)}^T\&times;\underset{\&OverBar;}{\mathrm{\&sigma;}},\underset{\&OverBar;}{\mathrm{\&epsiv;}}=\underset{\&OverBar;}{\mathrm{\&omega;}}\&times;(\underset{\&OverBar;}{J}\&CenterDot;\underset{\&OverBar;}{\mathrm{\&omega;}});$

步骤五、对变构型航天器的系统动力学方程积分，以变构型航天器所受外力作为己知约束条件，求得描述变构型航天器的动力学参量，建立变构型航天器的多刚体动力学模型。

所述的步骤一还包括以下子步骤：

步骤A、设变构型航天器最小可变结构单位为舱段，N记为航天器系统中舱段个数，每个舱段记作B_i(i＝1，…，N)，B₀表示变构型航天器以外运动为己知的假想舱段，固连在惯性坐标系中；相邻两个舱段间以铰连接，记作H_j，拓扑模型中采用一条连接邻接物体的有向线段表示；下标i、j分别表示物体和铰的标号；

采用规则标号法对舱段进行标号，规则标号法规定如下：

(1)与根物体B₀的邻接物体记为B₁，关联的铰为H₁；

(2)每个物体与其内接铰的序号相同；

(3)每个物体的序号大于其内接物体的序号；

(4)每个铰的指向一律背离B₀方向；

对于树系统，若B₀上有几组相互平行的树，即存在着相互独立的几个变构型航天器，则将每个变构型航天器看作一个独立的系统，分别进行建模；

对于非树系统，先按上述方法对其派生树系统进行标号，然后依次对切断铰补上标号，其方向由较小标号的物体指向较大标号的物体；

步骤B、根据变构型航天器的任务需求以及变型方案，将其变构型过程划分为n个变型阶段，每个变形阶段用邻接矩阵D_w(w＝1，…，n)表示航天器当前拓扑构型，定义n×n的邻接矩阵中的元素如下：

步骤C、根据邻接矩阵描述方法，取邻接矩阵D的上三角阵，记为则通路矩阵T满足ST＝1 _N，I _N表示N阶单位阵，1 _N为元素全为1的N阶列阵；

步骤D、定义内接物体矩阵L，用以描述变构型航天器拓扑结构信息，便于算法中的循环条件，提高运算效率，其形式如下：

L＝(L₁L₂…L_N)_R×N

其中L_i定义如下：关联矩阵S中第i行上所有非零元素的列标按降次排列，组成列向量L_i″，R＝max(L_i″行数)+1，维数不足R的L_i″在列尾补零，组成L_i。

内接物体矩阵L的行向量定义为R个N阶一维整型数组L^1-jj(i)(i＝1，…，N；jj＝1，…，R)，其中L¹(i)记为舱段B_i的内接物体数组，简写为L(i)。

所述的步骤一的步骤B中按照变胞方式的不同，可以通过以下三种形式得到D_w(w＝1，…，n)：

(1)舱段减少，其中D₀为初始构态，有n₀个舱段，变胞矩阵 为中间变量，采用二进制矩阵运算，即0+0＝0，1+0＝1，1+1＝0，具体表达形式如下：

$\begin{array}{cc}{U}_{(x,y),n}^p={\left(\begin{array}{cccccccc}1& 0& \mathrm{...}& 0& \mathrm{...}& 0& \mathrm{...}& 0\\ 0& 1& \mathrm{...}& 0& \mathrm{...}& 0& \mathrm{...}& 0\\ .& .& & .& & .& & .\\ .& .& & .& & .& & .\\ .& .& & .& & .& & .\\ 0& 0& \mathrm{...}& 1_{x,x}& \mathrm{...}& 1_{x,y}& \mathrm{...}& 0\\ .& .& & .& & .& & .\\ .& .& & .& & .& & .\\ .& .& & .& & .& & .\\ 0& 0& \mathrm{...}& 0& \mathrm{...}& 1_{y,y}& \mathrm{...}& 0\\ .& .& & .& & .& & .\\ .& .& & .& & .& & .\\ .& .& & .& & .& & .\\ 0& 0& \mathrm{...}& 0& \mathrm{...}& 0& \mathrm{...}& 1\end{array}\right)}_{n\&times;n}& E_{y,n}^p={\left(\begin{array}{ccccccc}1& 0& \mathrm{...}& 0& 0& \mathrm{...}& 0\\ 0& 1& \mathrm{...}& 0& 0& \mathrm{...}& 0\\ .& .& & .& .& & .\\ .& .& & .& .& & .\\ .& .& & .& .& & .\\ 0& 0& \mathrm{...}& 0_{y,y}& 0& \mathrm{...}& 0\\ 0& 0& \mathrm{...}& 0& 1& \mathrm{...}& 0\\ .& .& & .& .& & .\\ .& .& & .& .& & .\\ .& .& & .& .& & .\\ 0& 0& \mathrm{...}& 0& 0& \mathrm{...}& 1\end{array}\right)}_{n\&times;n}\end{array}$

(2)舱段增加， $\left\{\begin{array}{c}{D^{\&prime;}}_w=E_n{D}_0{E}_n^T\\ D_w={D^{\&prime;}}_w+\underset{p=w}{\overset{1}{\&Sigma;}}{W}_{(x,y),n_p}^p-\underset{p=w}{\overset{1}{\&Sigma;}}{V}_{(x,y),n_p}^p\end{array}\right.,$ 其中用来将原机构的邻接矩阵增加n-n₀行n-n₀列；表示增加的关联；表示断开的关联，形式同

$W_{(x,y),n_p}^p={\left(\begin{array}{ccccccc}0& \mathrm{...}& 0& \mathrm{...}& 0& \mathrm{...}& 0\\ .& & .& & .& & .\\ .& & .& & .& & .\\ .& & .& & .& & .\\ 0& \mathrm{...}& 1_{x,x}& \mathrm{...}& 0& \mathrm{...}& 0\\ .& & .& & .& & .\\ .& & .& & .& & .\\ .& & .& & .& & .\\ 0& \mathrm{...}& 0& \mathrm{...}& 1_{y,y}& \mathrm{...}& 0\\ .& & .& & .& & .\\ .& & .& & .& & .\\ .& & .& & .& & .\\ 0& \mathrm{...}& 0& \mathrm{...}& 0& \mathrm{...}& 0\end{array}\right)}_{n_p\&times;n_p}$

(3)舱段重组，任意构态p到p+1变化表示为式中意义同上，κ为系数，构件连接时κ＝1，断开时κ＝-1。

本发明的有益效果是：通过引入变胞理论描述方法，将传统建模前期需要人工进行的拓扑结构绘制准备变为可通过计算机自主进行的矩阵运算，弥补了变构型航天器对于地面实验及人工分析的过分依赖，增加了模型的通用性和程式化。本发明通过关联矩阵和内接物体矩阵的转化，弥补了变胞理论拓扑描述和运算难以建立数学运算关系，没有涉及系统的构态约束和功能需求，无法进一步准确描述系统的动力学特性以及不能与系统在特定环境下的应用紧密结合的缺陷。

附图说明

图1为本发明算法流程图；

图2为本发明邻接矩阵递推关系及参量示意图；

图3为具体实施案例的变构型航天器拓扑构型图；

图4为具体实施案例的连接铰广义角速度和角加速度变化曲线图。

具体实施方式

本发明包括以下步骤：

步骤一、根据变构型航天器的拓扑构型分析，建立变构型航天器的拓扑结构模型，获得变构型航天器的邻接矩阵D、关联矩阵S、内接物体矩阵L、内接物体数组L；

步骤二、建立坐标系，对变构型航天器进行运动学分析，包括以下步骤：

步骤A、建立参考惯性坐标系、体坐标系、铰坐标系，

建立右手惯性坐标系OXYZ，记作参考基e，基点O取在航天器初始构型D₀的整体质心处，X轴、Y轴、Z轴分别沿初始构型下的柯西惯性椭球长、中、短轴方向，满足右手坐标系定义；邻接两舱段中，舱段B_j为舱段B_i的内接舱段，即j＝L(i)；两舱段由铰H_i相连，在动力学模型中，忽略铰的形状和尺寸，简化为一个几何点，称为铰点，P与Q分别为铰点在舱段B_i与B_j上的位置；以舱段B_i与B_j的质心C_i与C_j为基点分别建立连体基，即舱段B_i与B_j的体坐标系，分别记作e ⁱ和e ^j，它们关于e的方向余弦阵分别记为A ⁱ与A ^j；建立以Q为原点的铰H_i的铰坐标系，即铰H_i的本地基，固结于点Q，记为过P建立铰H_i的动基e ^h，它固结于点P；

步骤B、采用广义坐标描述单舱段运动学参量，

定义广义坐标q _i(i＝1，…，Num)，组成列阵Num为铰的个数，根据规则标号法，对于树状拓扑结构的航天器系统，铰个数Num等于舱段个数N，刚体舱段B_i相对内接舱段B_j的相对位移相对速度转动角速度其中H^′hT、H^′ΩT是单位矢量动基e ^h的本地坐标阵，仅和铰的形式有关，统一用上标“′”表示相对于当地坐标系(即B_i浮动体坐标系)的矢量坐标阵，无上标表示相对于绝对坐标系的矢量坐标阵；

步骤三、根据单个变构型航天器舱段的运动学参量，建立邻接舱段运动学递推关系，描述变构型航天器运动，包括以下内容：

变构型航天器系统内各舱段角速度组成的列阵ω＝(ω₁，…，ω_N)^T，则有

$\underset{\&OverBar;}{\mathrm{\&omega;}}=\underset{\&OverBar;}{\mathrm{\&beta;}}\underset{\&OverBar;}{\stackrel{\&CenterDot;}{q}}+{\mathrm{\&omega;}}_0{\underset{\&OverBar;}{l}}_N,\underset{\&OverBar;}{\stackrel{\&CenterDot;}{\mathrm{\&omega;}}}=\underset{\&OverBar;}{\mathrm{\&beta;}}\underset{\&OverBar;}{\stackrel{\&CenterDot;\&CenterDot;}{q}}+\underset{\&OverBar;}{\mathrm{\&sigma;}}$

其中

各舱段质心位置列阵r＝(r₁，…，r_N)^T＝d ^*T 1 _N+(r₀+h₁)1 _N，质心速度其中h₁为六自由度虚铰H₁的滑移铰矢量，α＝-(H ^Ω T×d ^*+H ^h T)^T， $\underset{\&OverBar;}{\mathrm{\&upsi;}}=-{\underset{\&OverBar;}{d}}^{*T}\&times;{\mathrm{\&omega;}}_0{\underset{\&OverBar;}{l}}_N+({\mathrm{\&omega;}}_0\&times;h_1+{\stackrel{\&CenterDot;}{r}}_0){\underset{\&OverBar;}{l}}_N,$

d^*为广义通路矢量，其中c_ii为舱段B_i质心C_i指向B_i上铰点P的矢量；

步骤四、根据己知运动规律的舱段和步骤二中得到的邻接舱段递推关系，建立变构型航天器系统动力学方程，包括以下内容：

根据速度变分形式的动力学普遍方程，可得微分形式的系统动力学方程：

$\underset{\&OverBar;}{Z}\underset{\&OverBar;}{\stackrel{\&CenterDot;\&CenterDot;}{q}}=\underset{\&OverBar;}{z}$

其中Z为广义质量矩阵，z为广义力阵：

Z＝α ^*T mα ^*+β ^T Jβ

z＝α ^*T·(μF ^Out-mμw)+β ^T·(M ^O-J·σ-ε)+F ^eq+F ^τq

式中F ^Out为外力约束已知条件，F ^τq为直接作用于铰点的约束力元，F ^eq为作用于舱段上铰点以外的非约束力元，J为中心惯量张量J _i组成的对角阵，

$\underset{\&OverBar;}{\mathrm{\&mu;}}={\underset{\&OverBar;}{l}}_N-\frac{1}{m_s}\underset{\&OverBar;}{m}{\underset{\&OverBar;}{l}}_N{\underset{\&OverBar;}{l}}_N^T,{\underset{\&OverBar;}{\mathrm{\&alpha;}}}^{*}={\underset{\&OverBar;}{\mathrm{\&mu;}}}^T\underset{\&OverBar;}{\mathrm{\&alpha;}},{\underset{\&OverBar;}{\mathrm{\&mu;}}}^T\underset{\&OverBar;}{w}=-{\left({\underset{\&OverBar;}{d}}^{*}\underset{\&OverBar;}{\mathrm{\&mu;}}\right)}^T\&times;\underset{\&OverBar;}{\mathrm{\&sigma;}},\underset{\&OverBar;}{\mathrm{\&epsiv;}}=\underset{\&OverBar;}{\mathrm{\&omega;}}\&times;(\underset{\&OverBar;}{J}\&CenterDot;\underset{\&OverBar;}{\mathrm{\&omega;}});$

步骤五、对变构型航天器的系统动力学方程积分，以变构型航天器所受外力作为已知约束条件，求得描述变构型航天器的动力学参量，建立变构型航天器的多刚体动力学模型。

所述的步骤一还包括以下子步骤：

步骤A、设变构型航天器最小可变结构单位为舱段，N记为航天器系统中舱段个数，每个舱段记作B_O(i＝1，…，N)，B₀表示变构型航天器以外运动为已知的假想舱段，固连在惯性坐标系中：相邻两个舱段间以铰连接，记作H_j，拓扑模型中采用一条连接邻接物体的有向线段表示；下标i、j分别表示物体和铰的标号；

采用规则标号法对舱段进行标号，规则标号法规定如下：

(1)与根物体B₀的邻接物体记为B₁，关联的铰为H₁；

(2)每个物体与其内接铰的序号相同；

(3)每个物体的序号大于其内接物体的序号；

(4)每个铰的指向一律背离B₀方向；

对于树系统，若B₀上有几组相互平行的树，即存在着相互独立的几个变构型航天器，则将每个变构型航天器看作一个独立的系统，分别进行建模；

对于非树系统，先按上述方法对其派生树系统进行标号，然后依次对切断铰补上标号，其方向由较小标号的物体指向较大标号的物体；

步骤B、根据变构型航天器的任务需求以及变型方案，将其变构型过程划分为n个变型阶段，每个变形阶段用邻接矩阵D_w(w＝1，…，n)表示航天器当前拓扑构型，定义n×n的邻接矩阵中的元素如下：

步骤C、根据邻接矩阵描述方法，取邻接矩阵D的上三角阵，记为则通路矩阵T满足ST＝1 _N，I _N表示N阶单位阵，1 _N为元素全为1的N阶列阵；

步骤D、定义内接物体矩阵L，用以描述变构型航天器拓扑结构信息，便于算法中的循环条件，提高运算效率，其形式如下：

L＝(L₁L₂…L_N)_R×N

其中L_i定义如下：关联矩阵S中第i行上所有非零元素的列标按降次排列，组成列向量L_i″，R＝max(L_i″行数)+1，维数不足R的L_i″在列尾补零，组成L_i。

内接物体矩阵L的行向量定义为R个N阶一维整型数组L^1-jj(i)(i＝1，…，N；jj＝1，…，R)，其中L¹(i)记为舱段B_i的内接物体数组，简写为L(i)。

所述的步骤一的步骤B中按照变胞方式的不同，可以通过以下三种形式得到D_w(w＝1，…，n)：

(1)舱段减少，其中D₀为初始构态，有n₀个舱段，变胞矩阵 为中间变量，采用二进制矩阵运算，即0+0＝0，1+0＝1，1+1＝0，具体表达形式如下：

$\begin{array}{cc}{U}_{(x,y),n}^p={\left(\begin{array}{cccccccc}1& 0& \mathrm{...}& 0& \mathrm{...}& 0& \mathrm{...}& 0\\ 0& 1& \mathrm{...}& 0& \mathrm{...}& 0& \mathrm{...}& 0\\ .& .& & .& & .& & .\\ .& .& & .& & .& & .\\ .& .& & .& & .& & .\\ 0& 0& \mathrm{...}& 1_{x,x}& \mathrm{...}& 1_{x,y}& \mathrm{...}& 0\\ .& .& & .& & .& & .\\ .& .& & .& & .& & .\\ .& .& & .& & .& & .\\ 0& 0& \mathrm{...}& 0& \mathrm{...}& 1_{y,y}& \mathrm{...}& 0\\ .& .& & .& & .& & .\\ .& .& & .& & .& & .\\ .& .& & .& & .& & .\\ 0& 0& \mathrm{...}& 0& \mathrm{...}& 0& \mathrm{...}& 1\end{array}\right)}_{n\&times;n}& E_{y,n}^p={\left(\begin{array}{ccccccc}1& 0& \mathrm{...}& 0& 0& \mathrm{...}& 0\\ 0& 1& \mathrm{...}& 0& 0& \mathrm{...}& 0\\ .& .& & .& .& & .\\ .& .& & .& .& & .\\ .& .& & .& .& & .\\ 0& 0& \mathrm{...}& 0_{y,y}& 0& \mathrm{...}& 0\\ 0& 0& \mathrm{...}& 0& 1& \mathrm{...}& 0\\ .& .& & .& .& & .\\ .& .& & .& .& & .\\ .& .& & .& .& & .\\ 0& 0& \mathrm{...}& 0& 0& \mathrm{...}& 1\end{array}\right)}_{n\&times;n}\end{array}$

(2)舱段增加， $\left\{\begin{array}{c}{D^{\&prime;}}_w=E_n{D}_0{E}_n^T\\ D_w={D^{\&prime;}}_w+\underset{p=w}{\overset{1}{\&Sigma;}}{W}_{(x,y),n_p}^p-\underset{p=w}{\overset{1}{\&Sigma;}}{V}_{(x,y),n_p}^p\end{array}\right.,$ 其中用来将原机构的邻接矩阵增加n-n₀行n-n₀列；表示增加的关联；表示断开的关联，形式同

$W_{(x,y),n_p}^p={\left(\begin{array}{ccccccc}0& \mathrm{...}& 0& \mathrm{...}& 0& \mathrm{...}& 0\\ .& & .& & .& & .\\ .& & .& & .& & .\\ .& & .& & .& & .\\ 0& \mathrm{...}& 1_{x,x}& \mathrm{...}& 0& \mathrm{...}& 0\\ .& & .& & .& & .\\ .& & .& & .& & .\\ .& & .& & .& & .\\ 0& \mathrm{...}& 0& \mathrm{...}& 1_{y,y}& \mathrm{...}& 0\\ .& & .& & .& & .\\ .& & .& & .& & .\\ .& & .& & .& & .\\ 0& \mathrm{...}& 0& \mathrm{...}& 0& \mathrm{...}& 0\end{array}\right)}_{n_p\&times;n_p}$

(3)舱段重组，任意构态p到p+1变化表示为式中意义同上，κ为系数，构件连接时κ＝1，断开时κ＝-1。

下面结合附图和实施例对本发明进一步说明。

一种变构型航天器柔性多体动力学建模方法，其具体实施步骤包括：

首先，建立变构型航天器的拓扑结构模型。

如图3所示六板式变构型航天器，具有六个舱段，初始时所有连接铰松开，铰H₂在T_s时刻锁紧，当t≥T_s，锁紧效应等效为扭簧阻尼器作用。其余铰暂无驱动力矩，所有铰一旦到达目标角度，立刻锁紧，产生的锁紧力矩同铰H₂锁紧力矩。

采用规则标号法进行标号，舱段与舱段连接均用旋转铰，所有铰仅可围绕铰坐标系z轴(平行于体坐标系z轴)转动。

航天器变构型过程划分为2个变型阶段，初始阶段航天器为一个整体，即N＝1，邻接矩阵A₀＝[1]。之后航天器各舱段间电机工作，舱段数目N＝6，变构型航天器自由度改变。根据变胞理论，航天器邻接矩阵、关联矩阵和内接物体矩阵如下

根据得到的邻接舱段递推关系，建立变构型航天器系统动力学方程，然后将系统动力学方程降阶扩维成2N维一阶常微分方程组，在MATLAB中采用经典四阶龙格-库塔积分算法积分50s。

从附图4中可以看出，各舱段的运动在最初时刻存在较强的干扰，尤其是末端舱段的角速度变化很大，各舱段间出现较为严重的相互耦合。在时间到5秒左右，角速度曲线逐渐变为平稳，7秒到25秒为相对平滑段，这段各铰处在展开的中段。在第25秒到50秒积分结束，各铰角速度处在另一个相对平滑段，此段为各铰展开的末段，在空间表现为各构件相对平稳打开，趋于最终构型。

从各铰角度变化图上可以看出，0到5秒段角度变化不大，铰2和铰6相对有比较大的变化。但总体变化不大。这一时间段各舱段从折叠到打开，出现相互耦合现象，角度的变化有一定的反复。到8秒之后，角度变化趋势逐渐明显，表现为图上光滑变化曲线。可以看出，多刚体动力学建模比较复杂，在未考虑振动等耦合因素的情况下，连接铰的角速度变化在各舱段展开运动的初始阶段出现较强振荡。

上述分析证明了变构型航天器多刚体动力学建模方法能够表现出舱段位移速度等物理参量的传递效果，能够很好地描述出变构型航天器的动力学特性，并且可以给出航天器各个舱段实时姿态参数，为下一步进行航天器姿态控制的设计提供了较好的输入，使模型能够很好地与姿控系统接轨。

根据上述结果证明算法实现了利用变胞理论、邻接矩阵对多舱段航天器进行计算机自主拓扑的数学描述，并通过内接物体矩阵与动力学建模完好的结合，实现了程式化、通用性建模方法。

虽然结合了附图描述了本发明的实施方式，但对于本领域技术人员来说，在不脱离本发明原理的前提下，还可以做出若干改进，这些也应视为属于本发明的保护范围。

Claims (3)

1.一种变构型航天器多刚体动力学建模方法，其特征在于包括以下步骤：

步骤一、根据变构型航天器的拓扑构型分析，建立变构型航天器的拓扑结构模型，获得变构型航天器的邻接矩阵D、关联矩阵S、内接物体矩阵L、内接物体数组L；

步骤二、建立坐标系，对变构型航天器进行运动学分析，包括以下步骤：

步骤A、建立参考惯性坐标系、体坐标系、铰坐标系，

建立右手惯性坐标系OXYZ，记作参考基e，基点O取在航天器初始构型D₀的整体质心处，X轴、Y轴、Z轴分别沿初始构型下的柯西惯性椭球长、中、短轴方向，满足右手坐标系定义；邻接两舱段中，舱段B_j为舱段B_i的内接舱段，即j＝L(i)；两舱段由铰H_i相连，在动力学模型中简化为一个几何点，称为铰点，P与Q分别为铰点在舱段B_i与B_j上的位置；以舱段B_i与B_j的质心C_i与C_j为基点分别建立连体基，即舱段B_i与B_j的体坐标系，分别记作e ⁱ和e ^j，它们关于e的方向余弦阵分别记为A ⁱ与A ^j；建立以Q为原点的铰H_i的铰坐标系，即铰H_i的本地基，固结于点Q，记为过P建立铰H_i的动基e ^h，它固结于点P；

步骤B、采用广义坐标描述单舱段运动学参量，

定义广义坐标q _i(i＝1，…，Num)，组成列阵Num为铰的个数，根据规则标号法，对于树状拓扑结构的航天器系统，铰个数Num等于舱段个数N，刚体舱段B_i相对内接舱段B_j的相对位移相对速度转动角速度其中H ^′hT、H ^′ΩT是单位矢量动基e ^h的本地坐标阵，用上标表示相对于当地坐标系的矢量坐标阵，无上标表示相对于绝对坐标系的矢量坐标阵；

步骤三、根据单个变构型航天器舱段的运动学参量，建立邻接舱段运动学递推关系，描述变构型航天器运动，包括以下内容：

变构型航天器系统内各舱段角速度组成的列阵ω＝(ω₁，…，ω_N)^T，则有

$\underset{\&OverBar;}{\mathrm{\&omega;}}=\underset{\&OverBar;}{\mathrm{\&beta;}}\underset{\&OverBar;}{\stackrel{\&CenterDot;}{q}}\&times;{\mathrm{\&omega;}}_0{\underset{\&OverBar;}{l}}_N,\underset{\&OverBar;}{\stackrel{\&CenterDot;}{\mathrm{\&omega;}}}=\underset{\&OverBar;}{\mathrm{\&beta;}}\underset{\&OverBar;}{\stackrel{\&CenterDot;\&CenterDot;}{q}}+\underset{\&OverBar;}{\mathrm{\&sigma;}}$

其中β＝-TH ^ΩT，

各舱段质心位置列阵r＝(r₁，…，r_N)^T＝d ^*T l _N+(r₀+h₁)l _N，质心速度其中h₁为六自由度虚铰H₁的滑移铰矢量，α＝-(H ^Ω T×d ^*+H ^h T)^T， $\underset{\&OverBar;}{\mathrm{\&upsi;}}=-{\underset{\&OverBar;}{d}}^{*T}\&times;{\mathrm{\&omega;}}_0{\underset{\&OverBar;}{l}}_N+({\mathrm{\&omega;}}_0\&times;h_1+{\stackrel{\&CenterDot;}{r}}_0){\underset{\&OverBar;}{l}}_N,$

d^*为广义通路矢量，其中c_ii为舱段B_i质心C_i指向B_i上铰点P的矢量；

步骤四、根据已知运动规律的舱段和步骤二中得到的邻接舱段递推关系，建立变构型航天器系统动力学方程，包括以下内容：

根据速度变分形式的动力学普遍方程，可得微分形式的系统动力学方程：

$\underset{\&OverBar;}{Z}\underset{\&OverBar;}{\stackrel{\&CenterDot;\&CenterDot;}{q}}=\underset{\&OverBar;}{z}$

其中Z为广义质量矩阵，z为广义力阵：

Z＝α ^*T mα ^*+β ^T Jβ

z＝α ^*T·(μF ^Out-mμw)+β ^T·(M ^O-J·σ-ε)+F ^eq+F ^τq

式中F ^Out为外力约束已知条件，F ^τq为直接作用于铰点的约束力元，F ^eq为作用于舱段上铰点以外的非约束力元，J为中心惯量张量J _i组成的对角阵，

$\underset{\&OverBar;}{\mathrm{\&mu;}}={\underset{\&OverBar;}{l}}_N-\frac{1}{m_s}\underset{\&OverBar;}{m}{\underset{\&OverBar;}{l}}_N{\underset{\&OverBar;}{l}}_N^T,{\underset{\&OverBar;}{\mathrm{\&alpha;}}}^{*}={\underset{\&OverBar;}{\mathrm{\&mu;}}}^T\underset{\&OverBar;}{\mathrm{\&alpha;}},{\underset{\&OverBar;}{\mathrm{\&mu;}}}^T\underset{\&OverBar;}{w}=-{\left({\underset{\&OverBar;}{d}}^{*}\underset{\&OverBar;}{\mathrm{\&mu;}}\right)}^T\&times;\underset{\&OverBar;}{\mathrm{\&sigma;}},\underset{\&OverBar;}{\mathrm{\&epsiv;}}=\underset{\&OverBar;}{\mathrm{\&omega;}}\&times;(\underset{\&OverBar;}{J}\&CenterDot;\underset{\&OverBar;}{\mathrm{\&omega;}});$

步骤五、对变构型航天器的系统动力学方程积分，以变构型航天器所受外力作为已知约束条件，求得描述变构型航天器的动力学参量，建立变构型航天器的多刚体动力学模型。

2.根据权利要求1所述的变构型航天器多刚体动力学建模方法，其特征在于：所述的步骤一还包括以下子步骤：

步骤A、设变构型航天器最小可变结构单位为舱段，N记为航天器系统中舱段个数，每个舱段记作B_i(i＝1，…，N)，B₀表示变构型航天器以外运动为已知的假想舱段，固连在惯性坐标系中；相邻两个舱段间以铰连接，记作H_j，拓扑模型中采用一条连接邻接物体的有向线段表示；下标i、j分别表示物体和铰的标号；

采用规则标号法对舱段进行标号，规则标号法规定如下：

(1)与根物体B₀的邻接物体记为B₁，关联的铰为H₁；

(2)每个物体与其内接铰的序号相同；

(3)每个物体的序号大于其内接物体的序号；

(4)每个铰的指向一律背离B₀方向；

对于树系统，若B₀上有几组相互平行的树，即存在着相互独立的几个变构型航天器，则将每个变构型航天器看作一个独立的系统，分别进行建模；

对于非树系统，先按上述方法对其派生树系统进行标号，然后依次对切断铰补上标号，其方向由较小标号的物体指向较大标号的物体；

步骤B、根据变构型航天器的任务需求以及变型方案，将其变构型过程划分为n个变型阶段，每个变形阶段用邻接矩阵D_w(w＝1，…，n)表示航天器当前拓扑构型，定义n×n的邻接矩阵中的元素如下：

步骤C、根据邻接矩阵描述方法，取邻接矩阵D的上三角阵，记为则通路矩阵T满足ST＝l _N，l _N表示N阶单位阵，l _N为元素全为l的N阶列阵；

步骤D、定义内接物体矩阵L＝(L₁L₂…L_N)_R×N，其中L_i定义如下：关联矩阵S中第i行上所有非零元素的列标按降次排列，组成列向量L_i″，R＝max(L_i″行数)+1，维数不足R的L_i″在列尾补零，组成L_i；

内接物体矩阵L的行向量定义为R个N阶一维整型数组L^1-jj(i)(i＝1，…，N；jj＝1，…，R)，其中L^l(i)记为舱段B_i的内接物体数组，简写为L(i)。

3.根据权利要求2所述的变构型航天器多刚体动力学建模方法，其特征在于：所述的步骤一的步骤B中按照变胞方式的不同，可以通过以下三种形式得到D_w(w＝1，…，n)：

(1)舱段减少，其中D₀为初始构态，有n₀个舱段，变胞矩阵 为中间变量，采用二进制矩阵运算，即0+0＝0，1+0＝1，1+1＝0，具体表达形式如下：

$\begin{array}{cc}{U}_{(x,y),n}^p={\left(\begin{array}{cccccccc}1& 0& \mathrm{...}& 0& \mathrm{...}& 0& \mathrm{...}& 0\\ 0& 1& \mathrm{...}& 0& \mathrm{...}& 0& \mathrm{...}& 0\\ .& .& & .& & .& & .\\ .& .& & .& & .& & .\\ .& .& & .& & .& & .\\ 0& 0& \mathrm{...}& I_{x,x}& \mathrm{...}& I_{x,y}& \mathrm{...}& 0\\ .& .& & .& & .& & .\\ .& .& & .& & .& & .\\ .& .& & .& & .& & .\\ 0& 0& \mathrm{...}& 0& \mathrm{...}& I_{y,y}& \mathrm{...}& 0\\ .& .& & .& & .& & .\\ .& .& & .& & .& & .\\ .& .& & .& & .& & .\\ 0& 0& \mathrm{...}& 0& \mathrm{...}& 0& \mathrm{...}& 1\end{array}\right)}_{n\&times;n}& E_{y,n}^p={\left(\begin{array}{ccccccc}1& 0& \mathrm{...}& 0& 0& \mathrm{...}& 0\\ 0& 1& \mathrm{...}& 0& 0& \mathrm{...}& 0\\ .& .& & .& .& & .\\ .& .& & .& .& & .\\ .& .& & .& .& & .\\ 0& 0& \mathrm{...}& 0_{y,y}& 0& \mathrm{...}& 0\\ 0& 0& \mathrm{...}& 0& 1& \mathrm{...}& 0\\ .& .& & .& .& & .\\ .& .& & .& .& & .\\ .& .& & .& .& & .\\ 0& 0& \mathrm{...}& 0& 0& \mathrm{...}& 1\end{array}\right)}_{n\&times;n}\end{array}$

(2)舱段增加， $\left\{\begin{array}{c}{D^{\&prime;}}_w=E_n{D}_0{E}_n^T\\ D_w={D^{\&prime;}}_w+\underset{p=w}{\overset{1}{\&Sigma;}}{W}_{(x,y),n_p}^p-\underset{p=w}{\overset{1}{\&Sigma;}}{V}_{(x,y),n_p}^p\end{array}\right.,$ 其中 $E_n={\left[\begin{array}{cc}{\underset{\&OverBar;}{I}}_{n_0}& {\underset{\&OverBar;}{0}}_{n-n_0}\end{array}\right]}^T,$ 用来将原机构的邻接矩阵增加n-n₀行n-n₀列；表示增加的关联；表示断开的关联，形式同

$W_{(x,y),n_p}^p={\left(\begin{array}{ccccccc}0& \mathrm{...}& 0& \mathrm{...}& 0& \mathrm{...}& 0\\ .& & .& & .& & .\\ .& & .& & .& & .\\ .& & .& & .& & .\\ 0& \mathrm{...}& l_{x,x}& \mathrm{...}& 0& \mathrm{...}& 0\\ .& & .& & .& & .\\ .& & .& & .& & .\\ .& & .& & .& & .\\ 0& \mathrm{...}& 0& \mathrm{...}& l_{y,y}& \mathrm{...}& 0\\ .& & .& & .& & .\\ .& & .& & .& & .\\ .& & .& & .& & .\\ 0& \mathrm{...}& 0& \mathrm{...}& 0& \mathrm{...}& 0\end{array}\right)}_{n_p\&times;n_p}$

(3)舱段重组，任意构态p到p+1变化表示为式中意义同上，κ为系数，构件连接时κ＝1，断开时κ＝-1。