CN104808512A - 一种航天器多级驱动刚柔耦合响应的获取方法 - Google Patents

Abstract

一种航天器多级驱动刚柔耦合响应的获取方法，步骤为：(1)将航天器视为由本体、部件和关节三类结构组成；(2)分别建立本体、部件和关节的动能和势能；(3)根据连接条件，建立整个航天器的动能和势能；(4)利用步骤(3)得到的航天器动能和势能，采用Lagrange方程建立航天器的动力学方程；(5)求解步骤(4)得到的动力学方程，获得多级驱动过程中航天器的刚柔耦合动力学响应。本发明方法适用于具有柔性多体运动特征的复杂航天器，获取这类复杂多体结构航天器运动过程精确动力学响应，对控制系统设计和系统级仿真验证具有重要的意义。

Description

一种航天器多级驱动刚柔耦合响应的获取方法

技术领域

本发明属于航天器动力学领域，涉及一种针对柔性多体类航天器的动力学响应获取方法。

背景技术

随着载人航天和通信技术的发展，航天器正逐渐由单轴驱动太阳翼形式向着多体方向发展，例如载人空间机械臂、双轴驱动太阳翼、多自由度驱动天线等。这种刚体和柔性体的混合系统，在各体之间存在大范围的转动驱动运动的同时，还叠加有柔性体的弹性振动。这种复杂的刚柔耦合动力学效应，以及与控制系统的控制-结构交叉耦合，使得这类复杂航天器在动力学特性分析，特别是GNC控制方案验证时对柔性多体系统的动力学模型及其解算需求十分迫切。对于这种刚-柔混合多体系统，基于刚体假设的动力学模型已无法描述其复杂的动力学性态，必须采取柔性多体动力学建模方法，同时考虑部件大范围运动和部件本身的变形。

20世纪70年代中期，有关柔性多体系统动力学的理论研究工作实际上已经展开。到目前为止，柔性多体系统动力学的研究虽然取得了不少成果，但是还没有达到多刚体系统动力学的研究水平，其主要原因是对物体大范围运动与弹性变形耦合问题的认识与处理方法上遇到困难。而目前的商业软件都是封装形式，只能进行响应计算，而无法给出解析式的动力学方程及其中间计算过程，导致无法提供控制系统设计所需的解析式模型，以及无法开展动力学与控制的联合仿真。

发明内容

本发明解决的技术问题是：克服现有技术的不足，提供了一种航天器多级驱动刚柔耦合响应的获取方法，通过给出一套通用化的柔性多体动力学解析式模型，并对模型求解获取航天器的动力学响应，解决此类航天器复杂的机构运动与弹性振动的耦合动力学解算问题，可以为刚柔耦合多体结构航天器的动力学建模、刚柔耦合响应计算、控制系统设计，与控制系统的联合仿真提供技术支持。

本发明的技术解决方案是：一种航天器多级驱动刚柔耦合响应的获取方法，包括如下步骤：

(1)将航天器视为由本体、部件和关节三类结构组成，其中本体和部件之间，以及部件和部件之间通过关节连接，或者固定连接，部件为刚体或者柔性体，关节为带驱动功能的柔性体，包括电机和减速机构；

(2)获取本体的动能T_b，

$T_b=\frac{1}{2}{\left[\begin{array}{c}{v}_b\\ {\mathrm{\ω}}_b\end{array}\right]}^T{\stackrel{\‾}{M}}_b\left[\begin{array}{c}{v}_b\\ {\mathrm{\ω}}_b\end{array}\right],{\stackrel{\‾}{M}}_b=\left[\begin{array}{cc}{M}_b& P_b\\ P_b^T& I_b\end{array}\right]$

其中，v_b是本体参考系O_bX_bY_bZ_b相对于惯性坐标系O_iX_iY_iZ_i的速度，在O_bX_bY_bZ_b下描述，ω_b是本体参考系O_bX_bY_bZ_b相对于惯性坐标系O_iX_iY_iZ_i的角速度，在O_bX_bY_bZ_b下描述，M_b是本体的质量，P_b是本体相对于本体参考系O_bX_bY_bZ_b的静矩，I_b是本体相对于本体参考系O_bX_bY_bZ_b的转动惯量；

(3)获取部件的动能T_ai和势能V_ai，

$T_{\mathrm{ai}}=\frac{1}{2}{\left[\begin{array}{c}{\stackrel{\·}{\mathrm{\η}}}_{\mathrm{ai}}\\ {\stackrel{\·}{\mathrm{\ξ}}}_{\mathrm{ai},1}\\ {\stackrel{\·}{\mathrm{\ξ}}}_{\mathrm{ai},2}\end{array}\right]}^T{\stackrel{\‾}{M}}_{\mathrm{ai}}\left[\begin{array}{c}{\stackrel{\·}{\mathrm{\η}}}_{\mathrm{ai}}\\ {\stackrel{\·}{\mathrm{\ξ}}}_{\mathrm{ai},1}\\ {\stackrel{\·}{\mathrm{\ξ}}}_{\mathrm{ai},2}\end{array}\right],V_{\mathrm{ai}}=\frac{1}{2}{\left[\begin{array}{c}{\mathrm{\η}}_{\mathrm{ai}}\\ {\mathrm{\ξ}}_{\mathrm{ai},2}\end{array}\right]}^T{\stackrel{\‾}{K}}_{\mathrm{ai}}\left[\begin{array}{c}{\mathrm{\η}}_{\mathrm{ai}}\\ {\mathrm{\ξ}}_{\mathrm{ai},2}\end{array}\right]$

${\stackrel{\‾}{M}}_{\mathrm{ai}}=\left[\begin{array}{ccc}E& {\stackrel{\‾}{m}}_{\mathrm{ai},\mathrm{ij},1}+{\stackrel{\‾}{m}}_{\mathrm{ai},\mathrm{ij},2}T& {\stackrel{\‾}{m}}_{\mathrm{ai},\mathrm{ij},2}\\ {\stackrel{\‾}{m}}_{\mathrm{ai},\mathrm{ji},1}+T^T{\stackrel{\‾}{m}}_{\mathrm{ai},\mathrm{ji},2}& {\stackrel{\‾}{m}}_{\mathrm{ai},\mathrm{jj},11}+T^T{\stackrel{\‾}{m}}_{\mathrm{ai},\mathrm{jj},21}+{\stackrel{\‾}{m}}_{\mathrm{ai},\mathrm{jj},12}T+T^T{\stackrel{\‾}{m}}_{\mathrm{ai},\mathrm{jj},22}T& {\stackrel{\‾}{m}}_{\mathrm{ai},\mathrm{jj},12}+T^T{\stackrel{\‾}{m}}_{\mathrm{ai},\mathrm{jj},22}\\ {\stackrel{\‾}{m}}_{\mathrm{ai},\mathrm{ji},2}& {\stackrel{\‾}{m}}_{\mathrm{ai},\mathrm{jj},21}& {\stackrel{\‾}{m}}_{\mathrm{ai},\mathrm{jj},22}\end{array}\right]$

${\stackrel{\‾}{K}}_{\mathrm{ai}}=\left[\begin{array}{cc}{\stackrel{\‾}{k}}_{\mathrm{ai},\mathrm{ii}}& 0\\ 0& {\stackrel{\‾}{k}}_{\mathrm{ai},\mathrm{jj},22}\end{array}\right]$

其中，η_ai是部件i的正则模态集的模态坐标，ξ_ai,1是部件i的约束模态集相对于部件i-1连接位置的界面节点物理坐标，ξ_ai,2是部件i的约束模态集相对于部件i+1连接位置的界面节点物理坐标，和中的各矩阵是质量阵和刚度阵根据η_ai、ξ_ai,1和ξ_ai,2划分后的分块矩阵；

(4)获取关节的动能T_ri和势能V_ti，

$T_{\mathrm{ri}}=\frac{1}{2}{\left[\begin{array}{c}{v}_{\mathrm{ti}}\\ {\mathrm{\ω}}_{\mathrm{ti}}\end{array}\right]}^T{\stackrel{\‾}{M}}_{\mathrm{ii}}\left[\begin{array}{c}{v}_{\mathrm{ti}}\\ {\mathrm{\ω}}_{\mathrm{ti}}\end{array}\right],{\stackrel{\‾}{M}}_{\mathrm{ti}}=\left[\begin{array}{cc}{M}_{\mathrm{ti}}& P_{\mathrm{ti}}\\ P_{\mathrm{ti}}^T& I_{\mathrm{ti}}\end{array}\right]$

$V_{\mathrm{ti}}=\frac{1}{2}{\mathrm{\ξ}}_{\mathrm{ti}}^T{\stackrel{\‾}{K}}_{\mathrm{ti}}{\mathrm{\ξ}}_{\mathrm{ti}},{\stackrel{\‾}{K}}_{\mathrm{ti}}=\mathrm{diag}\{k_{\mathrm{ti},\mathrm{tx}},k_{\mathrm{ti},\mathrm{ty}},k_{\mathrm{ti},\mathrm{tz}},k_{\mathrm{ti},\mathrm{rx}},k_{\mathrm{ti},\mathrm{ry}},k_{\mathrm{ti},\mathrm{rz}}\}$

其中，v_ti是关节i电机转子参考系O_tiX_tiY_tiZ_ti相对于惯性系O_iX_iY_iZ_i的速度，ω_ti是关节i电机转子参考系O_tiX_tiY_tiZ_ti相对于惯性坐标系O_iX_iY_iZ_i的角速度，M_ti是关节i电机转子质量，P_ti是关节i电机转子相对于电机转子参考系O_tiX_tiY_tiZ_ti的静矩，I_ti是关节i电机转子相对于关节i电机转子参考系O_tiX_tiY_tiZ_ti的转动惯量，ξ_ti是关节i减速机构与电机转子连接点相对于关节i电机转子参考系O_tiX_tiY_tiZ_ti的弹性变形位移，K_ti是减速机构的刚度阵，k_ti,tx,k_ti,ty,k_ti,tz是3个平动刚度值，k_ti,rx,k_ti,ry,k_ti,rz是3个转动刚度值；

(5)根据连接条件，结合步骤(2)～(4)的结果，建立整个航天器的动能T和势能V，

$T={\mathrm{\ω}}^T\mathrm{I\ω}{\left[\begin{array}{c}{\mathrm{\ω}}_0\\ {\mathrm{\ω}}_1\\ \stackrel{\·}{\mathrm{\η}}\end{array}\right]}^T\left[\begin{array}{ccc}{I}_0\left(q\right)& I_{10}\left(q\right)& I_{20}\left(q\right)\\ I_{10}^T\left(q\right)& I_1\left(q\right)& I_{21}\left(q\right)\\ I_{20}^T\left(q\right)& I_{21}^T& I_2\left(q\right)\end{array}\right]\left[\begin{array}{c}{\mathrm{\ω}}_0\\ {\mathrm{\ω}}_1\\ \stackrel{\·}{\mathrm{\η}}\end{array}\right]$

$I_0=\left[\begin{array}{cc}{M}_s& P_s\\ P_s^T& I_s\end{array}\right],V={\mathrm{\η}}^T{\mathrm{\Ω}}^2\mathrm{\η},{\mathrm{\Ω}}^2=\mathrm{diag}({\mathrm{\Ω}}_1^2,{\mathrm{\Ω}}_2^2,...)$

其中M_s是航天器的质量，I_s是航天器相对于本体参考系O_bX_bY_bZ_b的转动惯量，P_s是航天器相对于本体参考系O_bX_bY_bZ_b的静矩， $\mathrm{\ω}={\left[\begin{array}{ccc}{\mathrm{\ω}}_0& {\mathrm{\ω}}_1& \stackrel{\·}{\mathrm{\η}}\end{array}\right]}^T,$ ω₀＝[v_b ω_b]^T是本体平动速度和转动角速度组成的向量， ${\mathrm{\ω}}_1={\left[\begin{array}{cccc}{\stackrel{\·}{q}}_1& {\stackrel{\·}{q}}_2& ...& {\stackrel{\·}{q}}_n\end{array}\right]}^T$ 是各部件或各关节的单轴转动角速度，q＝[q₁ q₂ ... q_n]^T是各部件或各关节的单轴转动角，为包含各部件的正则模态坐标、界面点坐标、关节弹性变形位移的速度项，I₀(q)、I₁(q)、I₂(q)、I₁₀(q)、I₂₀(q)、I₂₁(q)是广义质量阵I关于的分块矩阵；

所述的连接条件为：

$\left[\begin{array}{c}{v}_{t1}\\ {\mathrm{\ω}}_{t1}\end{array}\right]=\left[\begin{array}{c}{C}_{b\→t1}(v_b+{\tilde{r}}_{a1}^T{\mathrm{\ω}}_b)\\ {\stackrel{\‾}{\mathrm{\ω}}}_{t1}+C_{b\→t1}{\mathrm{\ω}}_b\end{array}\right]$

${\stackrel{\·}{\mathrm{\ξ}}}_{a\mathrm{1,1}}=\left[\begin{array}{c}{C}_{t1\→a1}(v_{t1}+{\tilde{r}}_{t1}^T{\mathrm{\ω}}_{t1})\\ C_{t1\→a1}{\mathrm{\ω}}_{t1}\end{array}\right]+C_{t1\→a1}{\stackrel{\·}{\mathrm{\ξ}}}_{t1}$

$\left[\begin{array}{c}{v}_{\mathrm{ti}}\\ {\mathrm{\ω}}_{\mathrm{ti}}\end{array}\right]=\left[\begin{array}{c}{C}_{\mathrm{ai}-1\→\mathrm{ti}}({\stackrel{\·}{\mathrm{\ξ}}}_{\mathrm{ai}-\mathrm{1,1}}^t+{\tilde{r}}_{\mathrm{ai}}^T{\stackrel{\·}{\mathrm{\ξ}}}_{\mathrm{ai}-\mathrm{1,1}}^r)\\ {\stackrel{\‾}{\mathrm{\ω}}}_{\mathrm{ti}}+C_{\mathrm{ai}-1\→\mathrm{ti}}{\stackrel{\·}{\mathrm{\ξ}}}_{\mathrm{ai}-\mathrm{1,1}}^r\end{array}\right]+C_{\mathrm{ai}-1\→\mathrm{ti}}{\stackrel{\·}{\mathrm{\ξ}}}_{\mathrm{ai}-\mathrm{1,2}},i>1$

${\stackrel{\·}{\mathrm{\ξ}}}_{\mathrm{ai},1}=\left[\begin{array}{c}{C}_{\mathrm{ti}\→\mathrm{ai}}(v_{\mathrm{ti}}+{\tilde{r}}_{\mathrm{ti}}^T{\mathrm{\ω}}_{\mathrm{ti}})\\ C_{\mathrm{ti}\→\mathrm{ai}}{\mathrm{\ω}}_{\mathrm{ti}}\end{array}\right]+C_{\mathrm{ti}\→\mathrm{ai}}{\stackrel{\·}{\mathrm{\ξ}}}_{\mathrm{ti}},i>1$

其中，r_a1是部件1与本体固定连接时部件1与本体连接点在O_bX_bY_bZ_b下的位置矢量，是r_a1的斜对称矩阵，是部件1与本体通过关节1连接时关节1参考系O_t1X_t1Y_t1Z_t1相对于本体参考系O_bX_bY_bZ_b的角速度，在O_t1X_t1Y_t1Z_t1下描述，r_t1是部件1与关节1连接点在O_t1X_t1Y_t1Z_t1下的位置矢量，是r_t1的斜对称矩阵，C_b→t1是O_bX_bY_bZ_b到O_t1X_t1Y_t1Z_t1的坐标变换矩阵，C_t1→a1是O_t1X_t1Y_t1Z_t1到部件1参考系O_a1X_a1Y_a1Z_a1的坐标变换矩阵，是关节i参考系O_tiX_tiY_tiZ_ti相对于部件i-1参考系O_ai-1X_ai-1Y_ai-1Z_ai-1的角速度，在O_tiX_tiY_tiZ_ti下描述，r_ai是部件i与部件i-1的连接点在O_ai-1X_ai-1Y_ai-1Z_ai-1下的位置矢量，是r_ai的斜对称矩阵，r_ti是部件i与关节i连接点在O_tiX_tiY_tiZ_ti下的位置矢量，是r_ti的斜对称矩阵，C_ai-1→ti是O_ai-1X_ai-1Y_ai-1Z_ai-1到O_tiX_tiY_tiZ_ti的坐标变换矩阵，C_ti→ai是O_tiX_tiY_tiZ_ti到O_aiX_aiY_aiZ_ai的坐标变换矩阵，是ξ_ai-1,1的平动分量，是ξ_ai-1,1的转动分量；

(6)利用步骤(5)得到的航天器的动能和势能，采用Lagrange方程建立航天器的动力学方程，

$\begin{array}{c}2I\stackrel{\·}{\mathrm{\ω}}+2\frac{\mathrm{dI}}{\mathrm{dt}}\mathrm{\ω}+2{\widetilde{\mathrm{\ω}}}_b\left[\begin{array}{ccc}{I}_0& I_{10}& I_{20}\\ 0& 0& 0\\ 0& 0& 0\end{array}\right]\mathrm{\ω}{\left[\begin{array}{ccc}0& {\mathrm{\ω}}^T\frac{\∂I}{\∂q_1}\mathrm{\ω}& {\mathrm{\ω}}^T\frac{\∂I}{\∂q_2}\mathrm{\ω}...0\end{array}\right]}^T+\\ \left[\begin{array}{c}0\\ 2{\tilde{v}}_b\\ 0\end{array}\left[\begin{array}{cc}E& 0\end{array}\right]{\tilde{I}}_0\mathrm{\ω}\right]+\left[\begin{array}{c}0\\ 0\\ 2\mathrm{\ξ\Ω}\stackrel{\·}{\mathrm{\η}}+{\mathrm{\Ω}}^2\mathrm{\η}\end{array}\right]=\left[\begin{array}{c}T\\ 0\end{array}\right]\end{array}$

$\frac{\mathrm{dI}}{\mathrm{dt}}=\frac{\∂I}{\∂q_1}{\stackrel{\·}{q}}_1+\frac{\∂I}{\∂q_2}{\stackrel{\·}{q}}_2+...\frac{\∂I}{\∂q_n}{\stackrel{\·}{q}}_n$

T＝[F_b T_b τ₁ τ₂ ... τ_n]^T

其中是ω_b的斜对称矩阵，是v_b的斜对称矩阵， ${\tilde{I}}_0=\left[\begin{array}{ccc}{I}_0& I_{10}& I_{20}\end{array}\right]$ 是广义质量阵I关于ω的分块矩阵，ξ是阻尼系数，Ω是振动频率矩阵，是I对时间t的导数，是I对q_i的偏导数，是q_i对时间t的导数，F_b,T_b是本体作用力和力矩，τ₁,τ₂,...τ_n是各部件电子转子驱动力矩；

(7)求解步骤(6)得到的动力学方程，获得航天器多级驱动过程中的刚柔耦合动力学响应。

本发明与现有技术相比的优点在于：本发明采用Lagrange方程的模态综合-混合坐标动力学建模方法。首先，将航天器划分为本体、部件和关节，分别建立其动能和势能；其次，根据本体、部件和关节之间的连接关系，建立整个航天器的动能和势能；最后，利用Lagrange方程，建立航天器动力学方程，积分求解获得航天器多级驱动刚柔耦合响应。本发明方法提供了用于控制系统设计的解析式动力学模型，并且能够精确的给出刚柔耦合响应结果，实现航天器多级驱动动力学与控制联合仿真。

附图说明

图1为本发明树形拓扑结构示意图；

图2为本发明柔性关节建模示意图；

图3为本发明方法的流程框图；

图4为本发明柔性体连接示意图。

具体实施方式

本发明所研究的多级驱动航天器是指航天器除本体外，带有多轴驱动附件，如多轴驱动天线、两轴太阳翼、空间机械臂等，这类航天器具有多体运动特征，同时部件和关节存在弹性变形。本文所指的关节为带驱动功能的机构，其主要由电机和减速机构组成。根据图1，本发明将航天器系统划分为本体和附件(如太阳翼、天线、机械臂等)；附件由若干个部件(如机械臂的臂杆)和关节组成，本体和部件，以及部件之间可以通过关节连接，或固定连接；部件可以考虑为刚体，或柔性体。本发明中，首先对部件、关节建模，其次给出本体、部件和关节之间的连接条件，再次推导出系统动能和势能，最后得到系统动力学方程并求解，得到航天器的耦合响应。

整个发明的实施流程如图3所示。具体实施方式详述如下：

(1)部件柔性建模

针对图3中的第i个部件，下面建立其动力学模型，其振动方程可以写为

$\stackrel{\·\·}{x}+\mathrm{kx}=f---\left(1\right)$

其中m为质量矩阵，k为刚度矩阵，f为包含有与第i-1个结构和第i+1个结构连接处的内力，x是节点响应。

当部件i与前后部件连接时，需要考虑其连接界面点，如图4中的ξ_ai,1和ξ_ai,2。因此，根据内部坐标和界面坐标，将方程(1)划分为

$\left[\begin{array}{cc}{m}_{\mathrm{ii}}& m_{\mathrm{ij}}\\ m_{\mathrm{ji}}& m_{\mathrm{jj}}\end{array}\right]\left[\begin{array}{c}{\stackrel{\·\·}{x}}_i\\ {\stackrel{\·\·}{x}}_j\end{array}\right]+\left[\begin{array}{cc}{k}_{\mathrm{ii}}& k_{\mathrm{ij}}\\ k_{\mathrm{ji}}& k_{\mathrm{jj}}\end{array}\right]\left[\begin{array}{c}{x}_i\\ x_j\end{array}\right]=\left[\begin{array}{c}0\\ f_j\end{array}\right]---\left(2\right)$

其中x_i为内部节点响应，x_j为界面节点响应，f_j为耦合界面连接内力，m_ii,m_ij,m_jj为质量矩阵m按照内、外部节点划分后的分块矩阵，k_ii,k_ij,k_jj为刚度矩阵k按照内、外部节点划分后的分块矩阵。

对于固定界面有x_j＝0，由方程(2)可获得子结构自由振动方程为

$m_{\mathrm{ii}}{\stackrel{\·\·}{x}}_i+k_{\mathrm{ii}}{x}_i=0---\left(3\right)$

于是，可求得子结构第n阶正则模态φ_n，满足

${\mathrm{\φ}}_n^T{m}_{\mathrm{ii}}{\mathrm{\φ}}_n=1---\left(4\right)$

${\mathrm{\φ}}_n^T{k}_{\mathrm{ii}}{\mathrm{\φ}}_n={\mathrm{\ω}}_n^2---\left(5\right)$

其中ω_n为固定界面第n阶频率。

假设φ_ii＝[φ₁,φ₂,…,φ_n]，则子结构的正则模态集为

${\mathrm{\Φ}}_n=\left[\begin{array}{c}{\mathrm{\φ}}_{\mathrm{ii}}\\ 0\end{array}\right]---\left(6\right)$

对式(2)，其静力方程为

$\left[\begin{array}{cc}{k}_{\mathrm{ii}}& k_{\mathrm{ij}}\\ k_{\mathrm{ji}}& k_{\mathrm{jj}}\end{array}\right]\left[\begin{array}{c}{x}_i\\ x_j\end{array}\right]=\left[\begin{array}{c}0\\ f_j\end{array}\right]---\left(7\right)$

由方程(7)可以得到

$x_i=-k_{\mathrm{ii}}^{-1}{k}_{\mathrm{ij}}{x}_j---\left(8\right)$

令得子结构的约束模态集

${\mathrm{\Φ}}_c=\left[\begin{array}{c}{\mathrm{\φ}}_{\mathrm{ij}}\\ E\end{array}\right]---\left(9\right)$

其中E是单位矩阵。

利用正则模态集Φ_n和约束模态集Φ_c组成子结构模态矩阵

Φ＝[Φ_n Φ_c]   (10)

将部件i物理坐标按照内外节点变换到模态坐标下有

x＝Φ_n,aiη_ai+Φ_c,ai，1ξ_ai，1+Φ_c,ai,2ξ_ai,2   (11)

其中，η_ai是部件i的正则模态集的模态坐标，ξ_ai,1是部件i的约束模态集相对于部件i-1连接位置的界面节点物理坐标，ξ_ai,2是部件i的约束模态集相对于部件i+1连接位置的界面节点物理坐标，Φ_n,ai是部件i的正则模态集，Φ_c,ai,1是部件i关于界面点ξ_ai,1的约束模态集，Φ_c,ai,2是部件i关于界面点ξ_ai,2的约束模态集。

那么，根据式(11)，将式(2)划分为

${\stackrel{\‾}{M}}_{\mathrm{ai}}\left[\begin{array}{c}{\stackrel{\·\·}{\mathrm{\η}}}_{\mathrm{ai}}\\ {\stackrel{\·\·}{\mathrm{\ξ}}}_{\mathrm{ai},1}\\ {\stackrel{\·\·}{\mathrm{\ξ}}}_{\mathrm{ai},2}\end{array}\right]+{\stackrel{\‾}{K}}_{\mathrm{ai}}\left[\begin{array}{c}{\mathrm{\η}}_{\mathrm{ai}}\\ {\mathrm{\ξ}}_{\mathrm{ai},1}\\ {\mathrm{\ξ}}_{\mathrm{ai},2}\end{array}\right]=0---\left(12\right)$

${\stackrel{\‾}{M}}_{\mathrm{ai}}=\left[\begin{array}{ccc}E& {\stackrel{\‾}{m}}_{\mathrm{ai},\mathrm{ij},1}& {\stackrel{\‾}{m}}_{\mathrm{ai},\mathrm{ij},2}\\ {\stackrel{\‾}{m}}_{\mathrm{ai},\mathrm{ji},1}& {\stackrel{\‾}{m}}_{\mathrm{ai},\mathrm{jj},11}& {\stackrel{\‾}{m}}_{\mathrm{ai},\mathrm{jj},12}\\ {\stackrel{\‾}{m}}_{\mathrm{ai},\mathrm{ji},2}& {\stackrel{\‾}{m}}_{\mathrm{ai},\mathrm{jj},21}& {\stackrel{\‾}{m}}_{\mathrm{ai},\mathrm{jj},22}\end{array}\right],{\stackrel{\‾}{K}}_{\mathrm{ai}}=\left[\begin{array}{ccc}{\stackrel{\‾}{k}}_{\mathrm{ai},\mathrm{ii}}& 0& 0\\ 0& {\stackrel{\‾}{k}}_{\mathrm{ai},\mathrm{jj},11}& {\stackrel{\‾}{k}}_{\mathrm{ai},\mathrm{jj},12}\\ 0& {\stackrel{\‾}{k}}_{\mathrm{ai},\mathrm{jj},21}& {\stackrel{\‾}{k}}_{\mathrm{ai},\mathrm{jj},22}\end{array}\right]---\left(13\right)$

其中和是模型(2)经过变换式(11)减缩后的质量阵和刚度阵，式(13)中的各矩阵是质量阵和刚度阵根据η_ai、ξ_ai,1和ξ_ai,2划分后的分块矩阵。

将界面点ξ_ai,2变换为与界面点ξ_ai,1的相对运动关系

ξ_ai,2＝Tξ_ai,1+ξ_ai,2   (14)

$T=\left[\begin{array}{cc}E& {-{\tilde{r}}_{\mathrm{ai}}}_{}\\ 0& E\end{array}\right]---\left(15\right)$

其中r_ai是部件i与部件i-1的连接位置矢量，是r_ai的斜对称矩阵。

那么，利用变换和分别左乘和右乘方程(12)

$\stackrel{\‾}{T}=\left[\begin{array}{ccc}E& 0& 0\\ 0& E& 0\\ 0& T& E\end{array}\right]---\left(16\right)$

变换后，最终获得柔性部件动力学方程(12)的系数矩阵

${\stackrel{\‾}{M}}_{\mathrm{ai}}=\left[\begin{array}{ccc}E& {\stackrel{\‾}{m}}_{\mathrm{ai},\mathrm{ij},1}+{\stackrel{\‾}{m}}_{\mathrm{ai},\mathrm{ij},2}T& {\stackrel{\‾}{m}}_{\mathrm{ai},\mathrm{ij},2}\\ {\stackrel{\‾}{m}}_{\mathrm{ai},\mathrm{ji},1}+T^T{\stackrel{\‾}{m}}_{\mathrm{ai},\mathrm{ji},2}& {\stackrel{\‾}{m}}_{\mathrm{ai},\mathrm{jj},11}+T^T{\stackrel{\‾}{m}}_{\mathrm{ai},\mathrm{jj},21}+{\stackrel{\‾}{m}}_{\mathrm{ai},\mathrm{jj},12}T+T^T{\stackrel{\‾}{m}}_{\mathrm{ai},\mathrm{jj},22}T& {\stackrel{\‾}{m}}_{\mathrm{ai},\mathrm{jj},12}+T^T{\stackrel{\‾}{m}}_{\mathrm{ai},\mathrm{jj},22}\\ {\stackrel{\‾}{m}}_{\mathrm{ai},\mathrm{ji},2}& {\stackrel{\‾}{m}}_{\mathrm{ai},\mathrm{jj},21}& {\stackrel{\‾}{m}}_{\mathrm{ai},\mathrm{jj},22}\end{array}\right]---\left(17\right)$

${\stackrel{\‾}{K}}_{\mathrm{ai}}=\left[\begin{array}{cc}{\stackrel{\‾}{k}}_{\mathrm{ai},\mathrm{ii}}& 0\\ 0& {\stackrel{\‾}{k}}_{\mathrm{ai},\mathrm{jj},22}\end{array}\right]---\left(18\right)$

当部件作为刚体考虑时，方程(12)退化为

$({\stackrel{\‾}{m}}_{\mathrm{ai},\mathrm{jj},11}+T^T{\stackrel{\‾}{m}}_{\mathrm{ai},\mathrm{jj},21}+{\stackrel{\‾}{m}}_{\mathrm{ai},\mathrm{jj},12}T+T^T{\stackrel{\‾}{m}}_{\mathrm{ai},\mathrm{jj},22}T){\stackrel{\·\·}{\mathrm{\ξ}}}_{\mathrm{ai},1}=0---\left(19\right)$

(2)关节柔性建模

对于由电机和减速机构组成的驱动机构(本发明中称为关节)，本发明将其模化为电机转子+弹簧系，其中电机转子为刚体，描述关节的质量特性，弹簧系描述减速机构的柔性，如图2所示，弹簧系由3个线弹簧和3个卷簧组成。在系统建模中，将电机转子按动能考虑为

$T_{\mathrm{ri}}=\frac{1}{2}{\left[\begin{array}{c}{v}_{\mathrm{ti}}\\ {\mathrm{\ω}}_{\mathrm{ti}}\end{array}\right]}^T{\stackrel{\‾}{M}}_{\mathrm{ti}}\left[\begin{array}{c}{v}_{\mathrm{ti}}\\ {\mathrm{\ω}}_{\mathrm{ti}}\end{array}\right],{\stackrel{\‾}{M}}_{\mathrm{ti}}=\left[\begin{array}{cc}{M}_{\mathrm{ti}}& P_{\mathrm{ti}}\\ P_{\mathrm{ti}}^T& I_{\mathrm{ti}}\end{array}\right]---\left(20\right)$

其中，v_ti是关节i电机转子参考系O_tiX_tiY_tiZ_ti相对于惯性系O_iX_iY_iZ_i的速度，ω_ti是电机转子参考系O_tiX_tiY_tiZ_ti相对于惯性坐标系O_iX_iY_iZ_i的角速度。M_ti是电机转子质量，P_ti是电机转子相对于其参考系O_tiX_tiY_tiZ_ti的静矩，I_ti是电机转子相对于其参考系O_tiX_tiY_tiZ_ti的转动惯量。电机转子参考系O_tiX_tiY_tiZ_ti的原点为电机转子连接端面的形心，三轴指向根据具体建模的需求进行定义。

在系统建模中，将电机转子按势能考虑为

$V_{\mathrm{ti}}=\frac{1}{2}{\mathrm{\ξ}}_{\mathrm{ti}}^T{\stackrel{\‾}{K}}_{\mathrm{ti}}{\mathrm{\ξ}}_{\mathrm{ti}},{\stackrel{\‾}{K}}_{\mathrm{ti}}=\mathrm{diag}\{k_{\mathrm{ti},\mathrm{tx}},k_{\mathrm{ti},\mathrm{ty}},k_{\mathrm{ti},\mathrm{tz}},k_{\mathrm{ti},\mathrm{rx}},k_{\mathrm{ti},\mathrm{ry}},k_{\mathrm{ti},\mathrm{rz}}\}---\left(21\right)$

其中，ξ_ti是图1中的弹簧系与驱动部件连接点相对于电机转子参考系的弹性变形位移，是弹簧系的刚度阵，k_ti,tx,k_ti,ty,k_ti,tz是3个平动刚度值，k_ti,rx,k_ti,ry,k_ti,rz是3个转动刚度值。

式(20)和式(21)为最终获得的柔性关节模型，当不考虑关节的柔性时，关节模型退化为式(20)。

(3)建立整个航天器的动能和势能

本发明所考虑的航天器力学模型如图1所示的树形拓扑结构。其中，O_iX_iY_iZ_i是惯性参考系，O_bX_bY_bZ_b是本体参考系，O_aiX_aiY_aiZ_ai是部件i参考系，O_tiX_tiY_tiZ_ti是关节i参考系。本体参考系O_bX_bY_bZ_b原点位于本体质心，三轴指向根据具体建模的需求进行定义。部件参考系O_aiX_aiY_aiZ_ai原点位于部件与前一个部件的连接点，三轴指向根据具体建模的需求进行定义。

采用Lagrange方程建立航天器动力学方程，下面首先建立各个体的动能和势能表达式，建模过程中不考虑重力影响。

本体动能表达式为

$T_b=\frac{1}{2}{\left[\begin{array}{c}{v}_b\\ {\mathrm{\ω}}_b\end{array}\right]}^T{\stackrel{\‾}{M}}_b\left[\begin{array}{c}{v}_b\\ {\mathrm{\ω}}_b\end{array}\right]---\left(22\right)$

${\stackrel{\‾}{M}}_b=\left[\begin{array}{cc}{M}_b& P_b\\ P_b^T& I_b\end{array}\right]---\left(23\right)$

其中，v_b是本体参考系O_bX_bY_bZ_b相对于惯性坐标系O_iX_iY_iZ_i的速度，在O_bX_bY_bZ_b下描述，ω_b是本体参考系O_bX_bY_bZ_b相对于惯性坐标系O_iX_iY_iZ_i的角速度，在O_bX_bY_bZ_b下描述。M_b是本体质量，P_b是本体相对于本体参考系O_bX_bY_bZ_b的静矩，I_b是本体相对于本体参考系O_bX_bY_bZ_b的转动惯量。

类似地，根据式(17)和式(18)，各部件动能和势能表达式为

$T_{\mathrm{ai}}=\frac{1}{2}{\left[\begin{array}{c}{\stackrel{\·}{\mathrm{\η}}}_{\mathrm{ai}}\\ {\stackrel{\·}{\mathrm{\ξ}}}_{\mathrm{ai},1}\\ {\stackrel{\·}{\mathrm{\ξ}}}_{\mathrm{ai},2}\end{array}\right]}^T{\stackrel{\‾}{M}}_{\mathrm{ai}}\left[\begin{array}{c}{\stackrel{\·}{\mathrm{\η}}}_{\mathrm{ai}}\\ {\stackrel{\·}{\mathrm{\ξ}}}_{\mathrm{ai},1}\\ {\stackrel{\·}{\mathrm{\ξ}}}_{\mathrm{ai},2}\end{array}\right]---\left(24\right)$

$V_{\mathrm{ai}}=\frac{1}{2}{\left[\begin{array}{c}{\mathrm{\η}}_{\mathrm{ai}}\\ {\mathrm{\ξ}}_{\mathrm{ai},2}\end{array}\right]}^T\left[\begin{array}{cc}{\stackrel{\‾}{k}}_{\mathrm{ai},\mathrm{ii}}& 0\\ 0& {\stackrel{\‾}{k}}_{\mathrm{ai},\mathrm{jj},22}\end{array}\right]\left[\begin{array}{c}{\mathrm{\η}}_{\mathrm{ai}}\\ {\mathrm{\ξ}}_{\mathrm{ai},2}\end{array}\right]---\left(25\right)$

其中界面点ξ_ai,1代表部件i的刚体运动，界面点ξ_ai,2代表部件i与部件i+1连接点相对于部件i参考系O_aiX_aiY_aiZ_ai的振动物理位移，代表部件i的正则模态坐标。

当考虑关节柔性时，其动能和势能表达式如式(20)和式(21)所示。

下面考虑本体、各部件和关节之间的连接条件。根据考虑关节柔性与否，本体和部件1的连接关系分为如下两种情况：

1)不考虑关节柔性

${\stackrel{\·}{\mathrm{\ξ}}}_{a1,1}=\left[\begin{array}{c}{C}_{b\→a1}(v_b+{\tilde{r}}_{a1}^T{\mathrm{\ω}}_b)\\ {\mathrm{\ω}}_{a1}{C}_{b\→a1}{\mathrm{\ω}}_b\end{array}\right]---\left(26\right)$

其中，ω_a1是部件1参考系O_a1X_a1Y_a1Z_a1相对于本体参考系O_bX_bY_bZ_b的角速度，在O_a1X_a1Y_a1Z_a1下描述，r_a1是部件1与本体连接点在O_bX_bY_bZ_b下的位置矢量，是r_a1的斜对称矩阵，C_b→a1是O_bX_bY_bZ_b到O_a1X_a1Y_a1Z_a1的坐标变换矩阵。

2)考虑关节柔性

$\left[\begin{array}{c}{v}_{t1}\\ {\mathrm{\ω}}_{t1}\end{array}\right]=\left[\begin{array}{c}{C}_{b\→t1}(v_b+{\tilde{r}}_{\mathrm{ai}}^T{\mathrm{\ω}}_b)\\ {\stackrel{\‾}{\mathrm{\ω}}}_{t1}+C_{b-\→t1}{\mathrm{\ω}}_b\end{array}\right]---\left(27\right)$

${\stackrel{\·}{\mathrm{\ξ}}}_{a1,1}=\left[\begin{array}{c}{C}_{t1\→a1}(v_{t1}+{\tilde{r}}_{t1}^T{\mathrm{\ω}}_{t1})\\ C_{t1\→a1}{\mathrm{\ω}}_{t1}\end{array}\right]+C_{t1\→a1}{\stackrel{\·}{\mathrm{\ξ}}}_{t1}---\left(28\right)$

其中，是关节1参考系O_t1X_t1Y_t1Z_t1相对于本体参考系O_bX_bY_bZ_b的角速度，在O_t1X_t1Y_t1Z_t1下描述，r_t1是部件1与关节1连接点在O_t1X_t1Y_t1Z_t1下的位置矢量，是r_t1的斜对称矩阵，C_b→t1是O_bX_bY_bZ_b到O_t1X_t1Y_t1Z_t1的坐标变换矩阵，C_t1→a1是O_t1X_t1Y_t1Z_t1到O_a1X_a1Y_a1Z_a1的坐标变换矩阵。

类似地，根据考虑关节柔性与否，部件i和部件i-1的连接关系分为如下两种情况：

1)不考虑关节柔性

${\stackrel{\·}{\mathrm{\ξ}}}_{\mathrm{ai},1}=\left[\begin{array}{c}{C}_{\mathrm{ai}-1\→\mathrm{ai}}({\stackrel{\·}{\mathrm{\ξ}}}_{\mathrm{ai}-\mathrm{1,1}}^t+{\tilde{r}}_{\mathrm{ai}}^T{\stackrel{\·}{\mathrm{\ξ}}}_{\mathrm{ai}-\mathrm{1,1}}^r)\\ {\stackrel{\‾}{\mathrm{\ω}}}_{\mathrm{ai}}+C_{\mathrm{ai}-1\→\mathrm{ai}}{\stackrel{\·}{\mathrm{\ξ}}}_{\mathrm{ai}-\mathrm{1,1}}^r\end{array}\right]+C_{\mathrm{ai}-1\→\mathrm{ai}}{\stackrel{\·}{\mathrm{\ξ}}}_{\mathrm{ai}-\mathrm{1,2}},i>1---\left(29\right)$

其中，是ξ_ai-1,1的平动分量，是ξ_ai-1,1的转动分量，ω_ai是部件i参考系O_aiX_aiY_aiZ_ai相对于部件i-1参考系O_ai-1X_ai-1Y_ai-1Z_ai-1的角速度,在O_aiX_aiY_aiZ_ai下描述，r_ai是部件i与部件i-1的连接点在O_ai-1X_ai-1Y_ai-1Z_ai-1下的位置矢量，是r_ai的斜对称矩阵，C_ai-1→ai是O_ai-1X_ai-1Y_ai-1Z_ai-1到O_aiX_aiY_aiZ_ai的坐标变换矩阵。

2)考虑关节柔性

$\left[\begin{array}{c}{v}_{\mathrm{ti}}\\ {\mathrm{\ω}}_{\mathrm{ti}}\end{array}\right]=\left[\begin{array}{c}{C}_{\mathrm{ai}-1\→\mathrm{ti}}({\stackrel{\·}{\mathrm{\ξ}}}_{\mathrm{ai}-\mathrm{1,1}}^t+{\tilde{r}}_{\mathrm{ai}}^T{\stackrel{\·}{\mathrm{\ξ}}}_{\mathrm{ai}-\mathrm{1,1}}^r)\\ {\stackrel{\‾}{\mathrm{\ω}}}_{\mathrm{ti}}+C_{\mathrm{ai}-1\→\mathrm{ti}}{\stackrel{\·}{\mathrm{\ξ}}}_{\mathrm{ai}-\mathrm{1,1}}^r\end{array}\right]+C_{\mathrm{ai}-1\→\mathrm{ti}}{\stackrel{\·}{\mathrm{\ξ}}}_{\mathrm{ai}-\mathrm{1,2}},i>1---\left(30\right)$

${\stackrel{\·}{\mathrm{\ξ}}}_{\mathrm{ai},1}=\left[\begin{array}{c}{C}_{\mathrm{ti}\→\mathrm{ai}}(v_{\mathrm{ti}}+{\tilde{r}}_{\mathrm{ti}}^T{\mathrm{\ω}}_{\mathrm{ti}})\\ C_{\mathrm{ti}\→\mathrm{ai}}{\mathrm{\ω}}_{\mathrm{ti}}\end{array}\right]+C_{\mathrm{ti}\→\mathrm{ai}}{\stackrel{\·}{\mathrm{\ξ}}}_{\mathrm{ti}},i>1---\left(31\right)$

其中，是关节i参考系O_tiX_tiY_tiZ_ti相对于部件i-1参考系O_ai-1X_ai-1Y_ai-1Z_ai-1的角速度，在O_tiX_tiY_tiZ_ti下描述，r_ti是部件i与关节i连接点在O_tiX_tiY_tiZ_ti下的位置矢量，是r_ti的斜对称矩阵，C_ai-1→ti是O_ai-1X_ai-1Y_ai-1Z_ai-1到O_tiX_tiY_tiZ_ti的坐标变换矩阵，C_ti→ai是O_tiX_tiY_tiZ_ti到O_aiX_aiY_aiZ_ai的坐标变换矩阵。

最后，根据连接条件(26)～(31)，得到系统的动能和势能表达式

T＝T_b+T_a1+T_a2+...+T_an  (32)

V＝V_a1+V_a2+...+V_an  (33)

(4)建立航天器动力学方程

下面根据(3)推导的系统动能方程和势能方程，利用Lagrange方程，推导系统动力学方程。将系统动能T和系统势能V描述如下

$T={\mathrm{\ω}}^T\mathrm{I\ω}{\left[\begin{array}{c}{\mathrm{\ω}}_0\\ {\mathrm{\ω}}_1\\ \stackrel{\·}{\mathrm{\η}}\end{array}\right]}^T\left[\begin{array}{ccc}{I}_0\left(q\right)& I_{10}\left(q\right)& I_{20}\left(q\right)\\ I_{10}^T\left(q\right)& I_1\left(q\right)& I_{21}\left(q\right)\\ I_{20}^T\left(q\right)& I_{21}^T& I_2\left(q\right)\end{array}\right]\left[\begin{array}{c}{\mathrm{\ω}}_0\\ {\mathrm{\ω}}_1\\ \stackrel{\·}{\mathrm{\η}}\end{array}\right]---\left(34\right)$

$V={\mathrm{\η}}^T{\mathrm{\Ω}}^2\mathrm{\η},{\mathrm{\Ω}}^2=\mathrm{diag}({\mathrm{\Ω}}_1^2,{\mathrm{\Ω}}_2^2,...)---\left(35\right)$

其中，系统坐标 $\mathrm{\ω}={\left[\begin{array}{ccc}{\mathrm{\ω}}_0& {\mathrm{\ω}}_1& \stackrel{\·}{\mathrm{\η}}\end{array}\right]}^T,$ ω₀＝[v_b ω_b]^T是本体平动速度和转动角速度， ${\mathrm{\ω}}_1={\left[\begin{array}{cccc}{\stackrel{\·}{q}}_1& {\stackrel{\·}{q}}_2& ...& {\stackrel{\·}{q}}_n\end{array}\right]}^T$ 是各部件或各关节的单轴转动角速度，q＝[q₁ q₂ ... q_n]^T是各部件或各关节的单轴转动角，为包含各部件的正则模态坐标、界面点坐标、关节弹性变形位移的速度项。Ω_i是第i阶振动频率矩阵。I₀(q)、I₁(q)、I₂(q)、I₁₀(q)、I₂₀(q)、I₂₁(q)是系统广义质量阵I关于的分块矩阵。

广义坐标和伪坐标形式的Lagrange方程如下

$\frac{d}{\mathrm{dt}}\frac{\∂L}{\∂{\stackrel{\·}{q}}_k}-\frac{\∂L}{\∂q_k}=Q_k,L=T-V,k=\mathrm{1,2},\·\·\·,n---\left(36\right)$

其中Q_k为对应全部主动力的广义力。

$(\frac{d}{\mathrm{dt}}+{\widetilde{\mathrm{\ω}}}_b)\frac{\∂T}{\∂{\mathrm{\ω}}_0}+\left[\begin{array}{c}0\\ {\tilde{v}}_b\frac{\∂T}{\∂v_b}\end{array}\right]=\left[\begin{array}{c}{F}_b\\ T_b\end{array}\right]---\left(37\right)$

其中，F_b为本体所受外力，T_b为本体所受外力矩，是ω_b的斜对称矩阵。

利用式(36)和式(37)，分别对式(34)和(35)中的本体伪坐标和广义坐标求解，得到系统动力学方程

$\begin{array}{c}2I\stackrel{\·}{\mathrm{\ω}}+2\frac{\mathrm{dI}}{\mathrm{dt}}\mathrm{\ω}+2{\widetilde{\mathrm{\ω}}}_b\left[\begin{array}{ccc}{I}_0& I_{10}& I_{20}\\ 0& 0& 0\\ 0& 0& 0\end{array}\right]\mathrm{\ω}{\left[\begin{array}{ccc}0& {\mathrm{\ω}}^T\frac{\∂I}{\∂q_1}\mathrm{\ω}& {\mathrm{\ω}}^T\frac{\∂I}{\∂q_2}\mathrm{\ω}...0\end{array}\right]}^T+\\ \left[\begin{array}{c}0\\ 2{\tilde{v}}_b\\ 0\end{array}\left[\begin{array}{cc}E& 0\end{array}\right]{\tilde{I}}_0\mathrm{\ω}\right]+\left[\begin{array}{c}0\\ 0\\ 2\mathrm{\ξ\Ω}\stackrel{\·}{\mathrm{\η}}+{\mathrm{\Ω}}^2\mathrm{\η}\end{array}\right]=\left[\begin{array}{c}T\\ 0\end{array}\right]\end{array}---\left(38\right)$

$\frac{\mathrm{dI}}{\mathrm{dt}}=\frac{\∂I}{\∂q_1}{\stackrel{\·}{q}}_1+\frac{\∂I}{\∂q_2}{\stackrel{\·}{q}}_2+...\frac{\∂I}{\∂q_n}{\stackrel{\·}{q}}_n---\left(39\right)$

T＝[F_b T_b τ₁ τ₂ ... τ_n]^T  (40)

其中是ω_b的斜对称矩阵，是v_b的斜对称矩阵， ${\tilde{I}}_0=\left[\begin{array}{ccc}{I}_0& I_{10}& I_{20}\end{array}\right],$ ξ是阻尼系数，Ω是振动频率矩阵，是I对时间t的导数，是I对q_i的偏导数，是q对时间t的导数，τ₁,τ₂,...τ_n是各部件驱动力矩。

综上所述，按照步骤(1)～(4)建立起来的方程(38)为最终获得的航天器多级驱动刚柔耦合动力学模型，对方程(38)进行数值积分运算，即可获得多级驱动过程中的刚柔耦合动力学响应。对方程(38)和控制系统进行编程实现，闭环求解即可实现航天器的动力学与控制联合仿真。

本发明说明书中未作详细描述的内容属本领域技术人员的公知技术。

Claims (1)

1.一种航天器多级驱动刚柔耦合响应的获取方法，其特征在于包括如下步骤：

(1)将航天器视为由本体、部件和关节三类结构组成，其中本体和部件之间，以及部件和部件之间通过关节连接，或者固定连接，部件为刚体或者柔性体，关节为带驱动功能的柔性体，包括电机和减速机构；

(2)获取本体的动能T_b，

$T_b=\frac{1}{2}{\left[\begin{array}{c}{v}_b\\ {\mathrm{\ω}}_b\end{array}\right]}^T{\stackrel{\‾}{M}}_b\left[\begin{array}{c}{v}_b\\ {\mathrm{\ω}}_b\end{array}\right],{\stackrel{\‾}{M}}_b=\left[\begin{array}{cc}{M}_b& P_b\\ {P_b}^T& I_b\end{array}\right]$

其中，v_b是本体参考系O_bX_bY_bZ_b相对于惯性坐标系O_iX_iY_iZ_i的速度，在O_bX_bY_bZ_b下描述，ω_b是本体参考系O_bX_bY_bZ_b相对于惯性坐标系O_iX_iY_iZ_i的角速度，在O_bX_bY_bZ_b下描述，M_b是本体的质量，P_b是本体相对于本体参考系O_bX_bY_bZ_b的静矩，I_b是本体相对于本体参考系O_bX_bY_bZ_b的转动惯量；

(3)获取部件的动能T_ai和势能V_ai，

$T_{\mathrm{ai}}=\frac{1}{2}{\left[\begin{array}{c}{\stackrel{.}{\mathrm{\η}}}_{\mathrm{ai}}\\ {\stackrel{.}{\mathrm{\ξ}}}_{\mathrm{ai},1}\\ {\stackrel{.}{\mathrm{\ξ}}}_{\mathrm{ai},2}\end{array}\right]}^T{\stackrel{\‾}{M}}_{\mathrm{ai}}\left[\begin{array}{c}{\stackrel{.}{\mathrm{\η}}}_{\mathrm{ai}}\\ {\stackrel{.}{\mathrm{\ξ}}}_{\mathrm{ai},1}\\ {\stackrel{.}{\mathrm{\ξ}}}_{\mathrm{ai},2}\end{array}\right],V_{\mathrm{ai}}=\frac{1}{2}{\left[\begin{array}{c}{\mathrm{\η}}_{\mathrm{ai}}\\ {\mathrm{\ξ}}_{\mathrm{ai},2}\end{array}\right]}^T{\stackrel{\‾}{K}}_{\mathrm{ai}}\left[\begin{array}{c}{\mathrm{\η}}_{\mathrm{ai}}\\ {\mathrm{\ξ}}_{\mathrm{ai},2}\end{array}\right]$

${\stackrel{\‾}{M}}_{\mathrm{ai}}=\left[\begin{array}{ccc}E& {\stackrel{\‾}{m}}_{\mathrm{ai},\mathrm{ij},1}+{\stackrel{\‾}{m}}_{\mathrm{ai},\mathrm{ij},2}T& {\stackrel{\‾}{m}}_{\mathrm{ai},\mathrm{ij},2}\\ {\stackrel{\‾}{m}}_{\mathrm{ai},\mathrm{ij},1}+T^T{\stackrel{\‾}{m}}_{\mathrm{ai},\mathrm{ji},2}& {\stackrel{\‾}{m}}_{\mathrm{ai},\mathrm{jj},11}+T^T{\stackrel{\‾}{m}}_{\mathrm{ai},\mathrm{jj},21}+{\stackrel{\‾}{m}}_{\mathrm{ai},\mathrm{jj},12}T+T^T{\stackrel{\‾}{m}}_{\mathrm{ai},\mathrm{jj},22}T& {\stackrel{\‾}{m}}_{\mathrm{ai},\mathrm{jj},12}+T^T{\stackrel{\‾}{m}}_{\mathrm{ai},\mathrm{jj},22}\\ {\stackrel{\‾}{m}}_{\mathrm{ai},\mathrm{ji},2}& {\stackrel{\‾}{m}}_{\mathrm{ai},\mathrm{jj},21}& {\stackrel{\‾}{m}}_{\mathrm{ai},\mathrm{jj},22}\end{array}\right]$

${\stackrel{\‾}{K}}_{\mathrm{ai}}=\left[\begin{array}{cc}{\stackrel{\‾}{k}}_{\mathrm{ai},\mathrm{ii}}& 0\\ 0& {\stackrel{\‾}{k}}_{\mathrm{ai},\mathrm{jj},22}\end{array}\right]$

其中，η_ai是部件i的正则模态集的模态坐标，ξ_ai,1是部件i的约束模态集相对于部件i-1连接位置的界面节点物理坐标，ξ_ai,2是部件i的约束模态集相对于部件i+1连接位置的界面节点物理坐标，和中的各矩阵是质量阵和刚度阵根据η_ai、ξ_ai,1和ξ_ai,2划分后的分块矩阵；

(4)获取关节的动能T_ri和势能V_ti，

$T_{\mathrm{ri}}=\frac{1}{2}{\left[\begin{array}{c}{v}_{\mathrm{ti}}\\ {\mathrm{\ω}}_{\mathrm{ti}}\end{array}\right]}^T{\stackrel{\‾}{M}}_{\mathrm{ti}}\left[\begin{array}{c}{v}_{\mathrm{ti}}\\ {\mathrm{\ω}}_{\mathrm{ti}}\end{array}\right],{\stackrel{\‾}{M}}_{\mathrm{ti}}=\left[\begin{array}{cc}{M}_{\mathrm{ti}}& P_{\mathrm{ti}}\\ P_{\mathrm{ti}}^T& I_{\mathrm{ti}}\end{array}\right]$

$V_{\mathrm{ti}}=\frac{1}{2}{\mathrm{\ξ}}_{\mathrm{ti}}^T{\stackrel{\‾}{K}}_{\mathrm{ti}}{\mathrm{\ξ}}_{\mathrm{yi}},{\stackrel{\‾}{K}}_{\mathrm{ti}}=\mathrm{diag}\{k_{\mathrm{ti},\mathrm{tx}},k_{\mathrm{ti},\mathrm{ty}},k_{\mathrm{ti},\mathrm{tz}},k_{\mathrm{ti},\mathrm{rx}},k_{\mathrm{ti},\mathrm{ry}},k_{\mathrm{ti},\mathrm{rz}}\}$

其中，v_ti是关节i电机转子参考系O_tiX_tiY_tiZ_ti相对于惯性系O_iX_iY_iZ_i的速度，ω_ti是关节i电机转子参考系O_tiX_tiY_tiZ_ti相对于惯性坐标系O_iX_iY_iZ_i的角速度，M_ti是关节i电机转子质量，P_ti是关节i电机转子相对于电机转子参考系O_tiX_tiY_tiZ_ti的静矩，I_ti是关节i电机转子相对于关节i电机转子参考系O_tiX_tiY_tiZ_ti的转动惯量，ξ_ti是关节i减速机构与电机转子连接点相对于关节i电机转子参考系O_tiX_tiY_tiZ_ti的弹性变形位移，是减速机构的刚度阵，k_ti,tx,k_ti,ty,k_ti,tz是3个平动刚度值，k_ti,rx,k_ti,ry,k_ti,rz是3个转动刚度值；

(5)根据连接条件，结合步骤(2)～(4)的结果，建立整个航天器的动能T和势能V，

$T={\mathrm{\ω}}^T\mathrm{I\ω}={\left[\begin{array}{c}{\mathrm{\ω}}_0\\ {\mathrm{\ω}}_1\\ \stackrel{.}{\mathrm{\η}}\end{array}\right]}^T\left[\begin{array}{ccc}{I}_0\left(q\right)& I_{10}\left(q\right)& I_{20}\left(q\right)\\ I_{10}^T\left(q\right)& I_1\left(q\right)& I_{21}\left(q\right)\\ I_{20}^T\left(q\right)& I_{21}^T\left(q\right)& I_2\left(q\right)\end{array}\right]\left[\begin{array}{c}{\mathrm{\ω}}_0\\ {\mathrm{\ω}}_1\\ \stackrel{.}{\mathrm{\η}}\end{array}\right]$

$I_0=\left[\begin{array}{cc}{M}_s& P_s\\ P_s^T& I_S\end{array}\right],v={\mathrm{\η}}^t{\mathrm{\Ω}}^2\mathrm{\η},{\mathrm{\Ω}}^2=\mathrm{diag}({\mathrm{\Ω}}_1^2{,\mathrm{\Ω}}_2^2,...)$

其中M_s是航天器的质量，I_s是航天器相对于本体参考系O_bX_bY_bZ_b的转动惯量，P_s是航天器相对于本体参考系O_bX_bY_bZ_b的静矩，ω₀＝[v_b ω_b]^T是本体平动速度和转动角速度组成的向量，是各部件或各关节的单轴转动角速度，q＝[q₁ q₂ ... q_n]^T是各部件或各关节的单轴转动角，为包含各部件的正则模态坐标、界面点坐标、关节弹性变形位移的速度项，I₀(q)、I₁(q)、I₂(q)、I₁₀(q)、I₂₀(q)、I₂₁(q)是广义质量阵I关于的分块矩阵；

所述的连接条件为：

$\left[\begin{array}{c}{v}_{t1}\\ {\mathrm{\ω}}_{t1}\end{array}\right]=\left[\begin{array}{c}{C}_{b\→t1}(v_b+{\tilde{r}}_{a1}^T{\mathrm{\ω}}_b)\\ {\stackrel{\‾}{\mathrm{\ω}}}_{t1}+C_{b\→t1}{\mathrm{\ω}}_b\end{array}\right]$

${\stackrel{.}{\mathrm{\ξ}}}_{a\mathrm{1,1}}=\left[\begin{array}{c}{C}_{t1\→a1}(v_{t1}+{\tilde{r}}_{t1}^T{\mathrm{\ω}}_{t1})\\ C_{t1\→a1}{\mathrm{\ω}}_{t1}\end{array}\right]+C_{t1\→a1}{\stackrel{.}{\mathrm{\ξ}}}_{t1}$

$\left[\begin{array}{c}{v}_{t1}\\ {\mathrm{\ω}}_{\mathrm{ti}}\end{array}\right]=\left[\begin{array}{c}{C}_{\mathrm{ai}\→1\→\mathrm{ti}}({\stackrel{.}{\mathrm{\ξ}}}_{\mathrm{ai}-\mathrm{1,1}}^t+{\tilde{r}}_{\mathrm{ai}}^T{\stackrel{.}{\mathrm{\ξ}}}_{\mathrm{ai}-\mathrm{1,1}}^r)\\ {\stackrel{\‾}{\mathrm{\ω}}}_{\mathrm{ti}}+C_{\mathrm{ai}-1\→\mathrm{ti}}{\stackrel{.}{\mathrm{\ξ}}}_{\mathrm{ai}-\mathrm{1,1}}^r\end{array}\right]+C_{\mathrm{ai}-1\→\mathrm{ti}}{\stackrel{.}{\mathrm{\ξ}}}_{\mathrm{ai}-\mathrm{1,2}},i>1$

${\stackrel{.}{\mathrm{\ξ}}}_{\mathrm{ai},1}=\left[\begin{array}{c}{C}_{\mathrm{ti}\→\mathrm{ai}}(v_{\mathrm{ti}}+{\tilde{r}}_{\mathrm{ti}}^T{\mathrm{\ω}}_{\mathrm{ti}})\\ C_{\mathrm{ti}\→\mathrm{ai}}{\mathrm{\ω}}_{\mathrm{ti}}\end{array}\right]+C_{\mathrm{ti}\→\mathrm{ai}}{\stackrel{.}{\mathrm{\ξ}}}_{\mathrm{ti}},i>1$

其中，r_a1是部件1与本体固定连接时部件1与本体连接点在O_bX_bY_bZ_b下的位置矢量，是r_a1的斜对称矩阵，是部件1与本体通过关节1连接时关节1参考系O_t1X_t1Y_t1Z_t1相对于本体参考系O_bX_bY_bZ_b的角速度，在O_t1X_t1Y_t1Z_t1下描述，r_t1是部件1与关节1连接点在O_t1X_t1Y_t1Z_t1下的位置矢量，是r_t1的斜对称矩阵，C_b→t1是O_bX_bY_bZ_b到O_t1X_t1Y_t1Z_t1的坐标变换矩阵，C_t1→a1是O_t1X_t1Y_t1Z_t1到部件1参考系O_a1X_a1Y_a1Z_a1的坐标变换矩阵，是关节i参考系O_tiX_tiY_tiZ_ti相对于部件i-1参考系O_ai-1X_ai-1Y_ai-1Z_ai-1的角速度，在O_tiX_tiY_tiZ_ti下描述，r_ai是部件i与部件i-1的连接点在O_ai-1X_ai-1Y_ai-1Z_ai-1下的位置矢量，是r_ai的斜对称矩阵，r_ti是部件i与关节i连接点在O_tiX_tiY_tiZ_ti下的位置矢量，是r_ti的斜对称矩阵，C_ai-1→ti是O_ai-1X_ai-1Y_ai-1Z_ai-1到O_tiX_tiY_tiZ_ti的坐标变换矩阵，C_ti→ai是O_tiX_tiY_tiZ_ti到O_aiX_aiY_aiZ_ai的坐标变换矩阵，是ξ_ai-1,1的平动分量，是ξ_ai-1,1的转动分量；

(6)利用步骤(5)得到的航天器的动能和势能，采用Lagrange方程建立航天器的动力学方程，

$\begin{array}{c}2I\stackrel{.}{\mathrm{\ω}}+2\frac{\mathrm{dI}}{\mathrm{dt}}\mathrm{\ω}+2{\widetilde{\mathrm{\ω}}}_b\left[\begin{array}{ccc}{I}_0& I_{10}& I_{20}\\ 0& 0& 0\\ 0& 0& 0\end{array}\right]\mathrm{\ω}-{\left[\begin{array}{ccc}0& {\mathrm{\ω}}^T\frac{\∂I}{{\∂q}_1}\mathrm{\ω}& {\mathrm{\ω}}^T\frac{\∂I}{{\∂q}_2}\mathrm{\ω}...0\end{array}\right]}^T+\\ \left[\begin{array}{c}0\\ 2{\tilde{v}}_b\left[\begin{array}{cc}E& 0\end{array}\right]\\ 0\end{array}\right]+\left[\begin{array}{c}0\\ 0\\ 2\mathrm{\ξ\Ω}\stackrel{.}{\mathrm{\η}}+{\mathrm{\Ω}}^2\mathrm{\η}\end{array}\right]=\left[\begin{array}{c}T\\ 0\end{array}\right]\end{array}$

$\frac{\mathrm{dI}}{\mathrm{dt}}=\frac{\∂I}{{\∂q}_1}{\stackrel{.}{q}}_1+\frac{\∂I}{{\∂q}_2}{\stackrel{.}{q}}_2+...\frac{\∂I}{{\∂q}_n}{\stackrel{.}{q}}_n$

T＝[F_b T_b τ₁ τ₂ ...τ_n]^T

其中是ω_b的斜对称矩阵，是v_b的斜对称矩阵，是广义质量阵I关于ω的分块矩阵，ξ是阻尼系数，Ω是振动频率矩阵，是I对时间t的导数，是I对q_i的偏导数，是q_i对时间t的导数，F_b,T_b是本体作用力和力矩，τ₁,τ₂,...τ_n是各部件电子转子驱动力矩；

(7)求解步骤(6)得到的动力学方程，获得航天器多级驱动过程中的刚柔耦合动力学响应。