CN102521440A - 一种基于积分插值法的空间绳系系统的仿真方法 - Google Patents

Abstract

本发明涉及一种基于积分插值法的空间绳系系统的仿真方法，采用系绳上点的应变ε来描述系绳的纵向运动，采用系绳上点的切向量τ来描述系绳的横向运动，并在模型离散化的过程中实现了ε和τ的分离求解，从而克服了系统的病态特性，提高了系统的求解效率；本发明在描述系绳的释放和回收过程时，也通过增加和减小第一段系绳的长度来说实现，但仍然按照柔性绳的特性来处理第一段系绳，避免了将系绳看作刚性杆所带来的误差，同时在插入新的节点时，本发明通过内插值的方法来计算新节点的状态，避免了主观设定新节点的状态所带来的仿真噪声。

Description

一种基于积分插值法的空间绳系系统的仿真方法

技术领域

本发明涉及一种基于积分插值法的空间绳系系统的仿真方法。

背景技术

空间绳系系统由于其作用范围大、机动灵活的特点而有着广泛的应用前景，目前，它已成为各航天大国空间技术发展的重点领域之一。在空间绳系系统的研究与应用过程中，对于系统数值仿真算法的研究是一项基础而非常重要的工作，但由于柔性系绳的存在，使得系统存在着明显的病态特性，同时系统在运动过程中系绳的长度、位置及形状的变化都会非常大，这与固体力学中所研究的刚体及挠性体有着很大不同，对于它们进行求解的方法也有着很大不同。

对于空间绳系系统的仿真，目前广泛使用的主要有3种算法：通过“珠子模型”直接构建系统运动的半离散方程，并利用“锚点”对系绳绳长和缩短的过程进行描述，然后通过相应的数值积分算法进行时间项的积分，这是一种直接差分的方法，算法简单，计算效率一般，计算精度也一般，另外“锚点”的引入也为系统带来了很大的仿真噪声；通过将系绳分成若干段，然后借助有限元的相关计算方法对连续形式的系统方程进行离散化求解。这种方法的精度非常高，能够对系绳的运动进行精细分析，但算法本身的复杂度非常高，计算效率低；通过选取满足系统几何边界条件的全局基函数，使用Galerkin或Ritz法对连续形式的系统方程进行离散求解。这种方法的计算效率高，计算速度快，但系统的计算精度依赖于主观选择的基函数。

在三种算法中，“珠子模型”由于其计算效率与精度均衡性较好，目前使用的非常多，有限元法较低的求解效率使得它近几年研究和使用的较少，而由于Galerkin或Ritz法的求解效率高，这种方法在近些年研究的较多，由此可以发现对于求解效率的重视是空间绳系系统仿真算法发展的趋势，因此如果能够给出一种方法，其计算精度与“珠子模型”相近，但计算效率却更高，这种算法不仅具有重要的学术价值，也具有广泛的应用前景。

发明内容

要解决的技术问题

为了避免现有技术的不足之处，本发明提出一种基于积分插值法的空间绳系系统的仿真方法，通过对于快速纵向运动分量和慢速横向运动分量的分离求解，使得在保持“珠子模型”精度的同时，提高了模型的求解效率。

技术方案

一种基于积分插值法的空间绳系系统的仿真方法，其特征在于步骤如下：

步骤1：建立空间绳系系统的运动模型

式中，s表示系绳上点的自然坐标，ξ表示释放点处的自然坐标，L表示系绳总的自然长度，ρ表示系绳密度，t表示系统时间，r表示空间中的点相对于母航天器B的矢径，n表示系绳上某点处的应力，f表示作用在系绳上的摄动力；

释放点的运动模型为：r(ξ，t)＝r_D(t)，

式中，r_D(t)表示关于时间的函数；

子航天器A的运动模型为：

式中，m_A表示子航天器A的质量，F表示作用在子航天器A上的摄动力；

绞盘式释放与回收机构的运动模型为：

式中，I₁和R₁分别表示绞盘机构的转动惯量及半径，

表示系绳绞盘上所缠绕的角度，M_C表示作用在绞盘上上的控制力矩；

所述

式中，μ_e表示地球引力常数，ω表示系统的轨道角速度，R_B表示母航天器B相对地球中心的矢径；

所述释放与回收过程中满足：

所述系绳上点的应变ε和系绳的切向量τ满足：(1+ε)τ＝r′；

所述系绳中应力n满足：

式中，E表示系绳的杨氏模量，A表示系绳的截面积，α表示系绳的粘性阻尼系数；

步骤2：使用N+1个节点将系绳分为N段，第一段系绳长度为η，以后各段系绳长度为h，并对所有节点依次进行编号，计算各个节点坐标的初值：

式中，

ε_i表示节点i的应变，τ_xi表示节点i处的切向量τ_i沿x轴方向的投影，τ_yi和τ_zi分别表示沿y和沿z轴的投影；

步骤3：构造系统状态微分

对于状态变量Z有：

得到函数f_z使得：

对于状态变量X有：

$\left\{\begin{array}{c}{u}_0={\left[(1+{\mathrm{\ϵ}}_0){\mathrm{\τ}}_0\right]}^{\·\·}=\frac{2}{\mathrm{\η}}g+\frac{2}{{\mathrm{\ρ\η}}^2}(n_1-n_0)+\frac{f_0+f_1}{\mathrm{\ρ\η}}-q_0\\ u_1={\left[(1+{\mathrm{\ϵ}}_1){\mathrm{\τ}}_1\right]}^{\·\·}=\frac{2}{\mathrm{\ρ\η}(h+\mathrm{\η})}(n_0-n_1)+\frac{2}{\mathrm{\ρh}(h+\mathrm{\η})}(n_2-n_1)+\frac{f_2-f_0}{\mathrm{\ρ}(h+\mathrm{\η})}-q_1\\ u_j={\left[(1+{\mathrm{\ϵ}}_j){\mathrm{\τ}}_j\right]}^{\·\·}=\frac{1}{{\mathrm{\ρh}}^2}(n_{j+1}-{2n}_j+n_{j-1})+\frac{f_{j+1}-f_{j-1}}{2\mathrm{\ρh}}-q_j(j=2,L,N-1)\\ u_N={\left[(1+{\mathrm{\ϵ}}_N){\mathrm{\τ}}_N\right]}^{\·\·}=-\frac{2}{{\mathrm{\ρh}}^2}(n_N-n_{N-1})-\frac{{2n}_N}{m_{A}h}-\frac{f_{N-1}+f_N}{\mathrm{\ρh}}+\frac{2F}{{\mathrm{hm}}_A}-q_N\end{array}\right.$

其中：

式中：

u_xj、u_yj和u_zj分别表示u_j沿x、y和z轴方向的投影，其中：

$A=\frac{1}{1+{\mathrm{\ϵ}}_j}\left[\begin{array}{cccc}{\mathrm{\τ}}_{\mathrm{yj}}^2+{\mathrm{\τ}}_{\mathrm{zj}}^2& -{\mathrm{\τ}}_{\mathrm{xj}}{\mathrm{\τ}}_{\mathrm{yj}}& -{\mathrm{\τ}}_{\mathrm{xj}}{\mathrm{\τ}}_{\mathrm{zj}}& (1+{\mathrm{\ϵ}}_j){\mathrm{\τ}}_{\mathrm{xj}}\\ -{\mathrm{\τ}}_{\mathrm{xj}}{\mathrm{\τ}}_{\mathrm{yj}}& {\mathrm{\τ}}_{\mathrm{xj}}^2+{\mathrm{\τ}}_{\mathrm{zj}}^2& -{\mathrm{\τ}}_{\mathrm{yj}}{\mathrm{\τ}}_{\mathrm{zj}}& (1+{\mathrm{\ϵ}}_j){\mathrm{\τ}}_{\mathrm{yj}}\\ -{\mathrm{\τ}}_{\mathrm{xj}}{\mathrm{\τ}}_{\mathrm{zj}}& -{\mathrm{\τ}}_{\mathrm{yj}}{\mathrm{\τ}}_{\mathrm{zj}}& {\mathrm{\τ}}_{\mathrm{xj}}^2+{\mathrm{\τ}}_{\mathrm{yj}}^2& (1+{\mathrm{\ϵ}}_j){\mathrm{\τ}}_{\mathrm{zj}}\\ (1+{\mathrm{\ϵ}}_j){\mathrm{\τ}}_{\mathrm{xj}}& (1+{\mathrm{\ϵ}}_j){\mathrm{\τ}}_{\mathrm{yj}}& (1+{\mathrm{\ϵ}}_j){\mathrm{\τ}}_{\mathrm{zj}}& -{(1+{\mathrm{\ϵ}}_j)}^2\end{array}\right]$

于是得到函数f_x和f_y使得：

将时间变量t_xz更新为t_xz+h_xz，当t_xz为h_y的整数倍时继续下一步，否则继续本步骤；

步骤4：使用四阶龙格-库塔方法求解X(t_xz+h_xz)和Z(t_xz+h_xz)：

$X(t_{\mathrm{xz}}+h_{\mathrm{xz}})=X\left(t_{\mathrm{xz}}\right)+\frac{h_{\mathrm{xz}}}{6}(k_{x1}+2k_{x2}+2k_{x3}+k_{x4})$

$Z(t_{\mathrm{xz}}+h_{\mathrm{xz}})=Z\left(t_{\mathrm{xz}}\right)+\frac{h_{\mathrm{xz}}}{6}(k_{z1}+2k_{z2}+2k_{z3}+k_{z4})$

其中：

k_x1＝h_xzf_x(t_xz，X(t_xz)，Y(t_y)，Z(t_xz))

k_z1＝h_xzf_z(t_xz，X(t_xz)，Y(t_y)，Z(t_xz))

$k_{x2}=h_{\mathrm{xz}}{f}_x(t_{\mathrm{xz}}+\frac{h_{\mathrm{xz}}}{2},X\left(t_{\mathrm{xz}}\right)+\frac{k_{x1}}{2},Y\left(t_y\right),Z\left(t_{\mathrm{xz}}\right)+\frac{k_{z1}}{2})$

$k_{z2}=h_{\mathrm{xz}}{f}_z(t_{\mathrm{xz}}+\frac{h_{\mathrm{xz}}}{2},X\left(t_{\mathrm{xz}}\right)+\frac{k_{x1}}{2},Y\left(t_y\right),Z\left(t_{\mathrm{xz}}\right)+\frac{k_{z1}}{2})$

$k_{x3}=h_{\mathrm{xz}}{f}_x(t_{\mathrm{xz}}+\frac{h_{\mathrm{xz}}}{2},X\left(t_{\mathrm{xz}}\right)+\frac{k_{x2}}{2},Y\left(t_y\right),Z\left(t_{\mathrm{xz}}\right)+\frac{k_{z2}}{2})$

$k_{z3}=h_{\mathrm{xz}}{f}_z(t_{\mathrm{xz}}+\frac{h_{\mathrm{xz}}}{2},X\left(t_{\mathrm{xz}}\right)+\frac{k_{x2}}{2},Y\left(t_y\right),Z\left(t_{\mathrm{xz}}\right)+\frac{k_{z2}}{2})$

k_x4＝h_xzf_x(t_xz+h_xz，X(t_xz)+k_x3，Y(t_y)，Z(t_xz)+k_z3)

k_z4＝h_xzf_z(t_xz+h_xz，X(t_xz)+k_x3，Y(t_y)，Z(t_xz)+k_z3)；

步骤5：就将时间变量t_xz更新为t_xz+h_xz，当t_xz为h_y的整数倍时进行下一步骤，否则重复步骤4；

步骤6：采用四阶龙格-库塔方法求解Y(t_y+h_xz)：

$Y(t_y+h_y)=Y\left(t_y\right)+\frac{h_y}{6}(k_{y1}+2k_{y2}+2k_{y3}+k_{y4})$

k_y1＝h_yf_y(t，X(t_y)，Y(t_y)，Z(t_y))

$k_{y2}=h_y{f}_y(t_y+\frac{h_y}{2},X(t_y+\frac{h_y}{2}),Y\left(t_y\right)+\frac{k_{y1}}{2},Z(t_{\mathrm{xz}}+\frac{h_y}{2}))$

$k_{y3}=h_y{f}_y(t_y+\frac{h_y}{2},X(t_y+\frac{h_y}{2}),Y\left(t_y\right)+\frac{k_{y2}}{2},Z(t_y+\frac{h_y}{2}))$

k_y4＝h_yf_y(t_y+h_y，X(t_y+h_y)，Y(t_y)+k_y3，Z(t_y+h_y))

当系绳分段数N达到分段上限SupLim时，将柔性绳上编号为奇数的节点删除，使相邻分段合并；所述SupLim为偶数，范围为

当系绳分段数N达到分段下限InfLim时，则在节点i和节点i+1(i＝1，2，K，N-1)的中点处插入新的节点

其节点状态由节点节点i和节点i+1的状态取平均得到；所述InfLim范围为[1，100]；

当η≥1.5h，在节点0和节点1之间插入一个新的节点，使得新的节点与节点1之间绳段的自然长度为h，且新节点的状态由线性插值获得：

${\mathrm{\ϵ}}_{*}=\frac{h}{\mathrm{\η}}{\mathrm{\ϵ}}_0+\frac{\mathrm{\η}-h}{\mathrm{\η}}{\mathrm{\ϵ}}_1,$ ${\mathrm{\τ}}_{*}=\frac{h}{\mathrm{\η}}{\mathrm{\τ}}_0+\frac{\mathrm{\η}-h}{\mathrm{\η}}{\mathrm{\τ}}_1,$

当η＜0.5h，将第一段与第二段系绳合并起来且删去节点1；

步骤7：将时间变量t_y更新为t_y+h_y，若t_y达到仿真结束的时间则结束仿真；若未达到，则从步骤3重新开始。

有益效果

本发明提出的一种基于积分插值法的空间绳系系统的仿真方法，与现有技术相比，本发明的有益效果是：

一、柔性系绳中的运动可以分解为纵向运动和横向运动，而纵向运动比横向运动的传播速度快得多，这使得空间绳系系统有着明显的病态特性。本发明在计算过程中，采用系绳上点的应变ε来描述系绳的纵向运动，采用系绳上点的切向量τ来描述系绳的横向运动，并在模型离散化的过程中实现了ε和τ的分离求解，从而克服了系统的病态特性，提高了系统的求解效率；

二、“珠子模型”将释放点与第一个“珠点”之间的系绳假设为刚性杆，通过增加和减小刚性杆的长度来实现系绳的释放和回收，当刚性杆的长度增加或减小到一定程度时，就在刚性杆上插入一个新的“珠点”或者删去第一个“珠点”，而新的“珠点”的状态是主观设定的，这给“珠子模型”的计算带来了很大的噪声。本发明在描述系绳的释放和回收过程时，也通过增加和减小第一段系绳的长度来说实现，但仍然按照柔性绳的特性来处理第一段系绳，避免了将系绳看作刚性杆所带来的误差，同时在插入新的节点时，本发明通过内插值的方法来计算新节点的状态，避免了主观设定新节点的状态所带来的仿真噪声。

附图说明

图1：为空间绳系系统的示意图；

图2：为绞盘式释放与回收机构的示意图；

图3：为本发明采用的分段方式示意图；

图4：为本发明采用的分段合并方式示意图；

图5：为本发明采用的分段分裂方式示意图；

图6：为本发明采用的节点插入方式示意图；

图7：为本发明采用的节点删除方式示意图；

图8：为本发明总体思路的流程图；

具体实施方式

现结合实施例、附图对本发明作进一步描述：

(1)建立系统数学模型：

对于如图1所示的空间绳系系统，系绳的运动可以描述为：

式中，系绳总的自然长度L＝500m，系绳的线密度ρ＝4.524×10^-3kg/m，系绳上的摄动力f＝0，且

式中，地球引力常数μ_e＝398600.5km³/s²，母航天器B的轨道为标准圆形，轨道半径R_B＝6871km，R_B＝(6.871×10⁶，0，0)^T，ω_B＝0.0011085rad/s。

对于释放点，它的运动可以描述为：

r(ξ，t)＝0

对于子航天器A，它的运动可以描述为：

式中，子航天器的质量m_A＝10kg，作用在子航天器上的摄动力F＝0。

采用恒速释放，故：

另外，释放点处的自然坐标满足：

式中，R₁＝0.2m。系绳上点的应变ε和系绳的切向量τ满足：

(1+ε)τ＝r′

系绳中应力n满足：

式中，系绳的弹性刚度EA＝1.04×10⁵N，系绳的粘性阻尼系绳α＝0。

(2)系绳初始处于沿径向的无形变的伸直状态，初始长度为4m，如图所示。初始将系绳均匀分为4段，第一段长度η＝1m，以后三段长度h＝1m，根据系统的初始状态，可以计算系统初值：

X(0)＝(0，0，0，0，L，0，0)

Y(0)＝(1，0，0，0，0，0，L，1，0，0，0，0，0)

Z(0)＝(2480，5，496)

另外，在计算开始前将时间变量t_xz和t_y都置为0。

(3)构造系统状态微分：

对于状态变量Z有：

于是可以得到函数f_z使得：

对于状态变量X有：

$\left\{\begin{array}{c}{u}_0={\left[(1+{\mathrm{\ϵ}}_0){\mathrm{\τ}}_0\right]}^{\·\·}=\frac{2}{\mathrm{\η}}g+\frac{2}{{\mathrm{\ρ\η}}^2}(n_1-n_0)+\frac{f_0+f_1}{\mathrm{\ρ\η}}-q_0\\ u_1={\left[(1+{\mathrm{\ϵ}}_1){\mathrm{\τ}}_1\right]}^{\·\·}=\frac{2}{\mathrm{\ρ\η}(h+\mathrm{\η})}(n_0-n_1)+\frac{2}{\mathrm{\ρh}(h+\mathrm{\η})}(n_2-n_1)+\frac{f_2-f_0}{\mathrm{\ρ}(h+\mathrm{\η})}-q_1\\ u_j={\left[(1+{\mathrm{\ϵ}}_j){\mathrm{\τ}}_j\right]}^{\·\·}=\frac{1}{{\mathrm{\ρh}}^2}(n_{j+1}-{2n}_j+n_{j-1})+\frac{f_{j+1}-f_{j-1}}{2\mathrm{\ρh}}-q_j(j=2,L,N-1)\\ u_N={\left[(1+{\mathrm{\ϵ}}_N){\mathrm{\τ}}_N\right]}^{\·\·}=-\frac{2}{{\mathrm{\ρh}}^2}(n_N-n_{N-1})-\frac{{2n}_N}{m_{A}h}-\frac{f_{N-1}+f_N}{\mathrm{\ρh}}+\frac{2F}{{\mathrm{hm}}_A}-q_N\end{array}\right.$ 式中，

另外有

式中，u_xj、u_yj和u_zj分别表示u_j沿x、y和z轴方向的投影，

$A=\frac{1}{1+{\mathrm{\ϵ}}_j}\left[\begin{array}{cccc}{\mathrm{\τ}}_{\mathrm{yj}}^2+{\mathrm{\τ}}_{\mathrm{zj}}^2& -{\mathrm{\τ}}_{\mathrm{xj}}{\mathrm{\τ}}_{\mathrm{yj}}& -{\mathrm{\τ}}_{\mathrm{xj}}{\mathrm{\τ}}_{\mathrm{zj}}& (1+{\mathrm{\ϵ}}_j){\mathrm{\τ}}_{\mathrm{xj}}\\ -{\mathrm{\τ}}_{\mathrm{xj}}{\mathrm{\τ}}_{\mathrm{yj}}& {\mathrm{\τ}}_{\mathrm{xj}}^2+{\mathrm{\τ}}_{\mathrm{zj}}^2& -{\mathrm{\τ}}_{\mathrm{yj}}{\mathrm{\τ}}_{\mathrm{zj}}& (1+{\mathrm{\ϵ}}_j){\mathrm{\τ}}_{\mathrm{yj}}\\ -{\mathrm{\τ}}_{\mathrm{xj}}{\mathrm{\τ}}_{\mathrm{zj}}& -{\mathrm{\τ}}_{\mathrm{yj}}{\mathrm{\τ}}_{\mathrm{zj}}& {\mathrm{\τ}}_{\mathrm{xj}}^2+{\mathrm{\τ}}_{\mathrm{yj}}^2& (1+{\mathrm{\ϵ}}_j){\mathrm{\τ}}_{\mathrm{zj}}\\ (1+{\mathrm{\ϵ}}_j){\mathrm{\τ}}_{\mathrm{xj}}& (1+{\mathrm{\ϵ}}_j){\mathrm{\τ}}_{\mathrm{yj}}& (1+{\mathrm{\ϵ}}_j){\mathrm{\τ}}_{\mathrm{zj}}& -{(1+{\mathrm{\ϵ}}_j)}^2\end{array}\right]$

于是可以得到函数f_x和f_y使得：

(4)对于状态变量X和Z，选取步长h_xz＝1×10⁴s，使用四阶龙格-库塔方法计算下一步的状态X(t_xz+h_xz)和Z(t_xz+h_xz)，在计算过程中认为Y保持在Y(t_y)。

(5)将时间变量t_xz更新为t_xz+h_xz，选取状态变量Y的计算步长h_y＝0.001s，并判断t_xz是否为h_y的整数倍，如果是，则转至第6步，如果不是，则转至第4步继续进行计算。

(6)使用四阶龙格-库塔方法计算状态变量Y下一步的状态Y(t_y+h_xz)：

$Y(t_y+h_y)=Y\left(t_y\right)+\frac{h_y}{6}(k_{y1}+2k_{y2}+2k_{y3}+k_{y4})$

k_y1＝h_yf_y(t，X(t_y)，Y(t_y)，Z(t_y))

$k_{y2}=h_y{f}_y(t_y+\frac{h_y}{2},X(t_y+\frac{h_y}{2}),Y\left(t_y\right)+\frac{k_{y1}}{2},Z(t_{\mathrm{xz}}+\frac{h_y}{2}))$

$k_{y3}=h_y{f}_y(t_y+\frac{h_y}{2},X(t_y+\frac{h_y}{2}),Y\left(t_y\right)+\frac{k_{y2}}{2},Z(t_y+\frac{h_y}{2}))$

k_y4＝h_yf_y(t_y+h_y，X(t_y+h_y)，Y(t_y)+k_y3，Z(t_y+h_y))

若系绳的分段数N达到分段上限SupLim＝150，则将柔性绳上编号为奇数的节点删除，使相邻分段合并，如图4所示。

若系绳的分段数N达到分段下限InfLim＝1(释放过程，下限直接设为1)，则在节点i和节点i+1(i＝1，2，K，N-1)的中点处插入新的节点

其节点状态由节点节点i和节点i+1的状态取平均得到，如图5所示。

若η≥1.5h，则在节点0和节点1之间插入一个新的节点，使得新的节点与节点1之间绳段的自然长度为h，如图6所示，且新节点的状态可由线性插值获得，即

${\mathrm{\ϵ}}_{*}=\frac{h}{\mathrm{\η}}{\mathrm{\ϵ}}_0+\frac{\mathrm{\η}-h}{\mathrm{\η}}{\mathrm{\ϵ}}_1,$ ${\mathrm{\τ}}_{*}=\frac{h}{\mathrm{\η}}{\mathrm{\τ}}_0+\frac{\mathrm{\η}-h}{\mathrm{\η}}{\mathrm{\τ}}_1,$

若η＜0.5h，则将第一段与第二段系绳合并起来，删去节点1，如图7所示。

(7)将时间变量t_y更新为t_y+h_y，同时判断t_y是否达到300s，若达到则结束仿真，若未达到，则转至第3步。

仿真试验

本发明的效果可以通过以下仿真实验结果进行验证：

在同一台主频为2.8GHz计算机(CPU为Core E5500)上，同时运行由VC2010所编写的本发明方法和经典“珠子模型”的求解程序，表1中给出了两种方法求解结果的对比。见表1

由表1可以看出：在基本保持“珠子模型”计算精度的同时，积分插值法实现了计算效率的大幅度提升。

Claims (1)

1.一种基于积分插值法的空间绳系系统的仿真方法，其特征在于步骤如下：

步骤1：建立空间绳系系统的运动模型

式中，s表示系绳上点的自然坐标，ξ表示释放点处的自然坐标，L表示系绳总的自然长度，ρ表示系绳密度，t表示系统时间，r表示空间中的点相对于母航天器B的矢径，n表示系绳上某点处的应力，f表示作用在系绳上的摄动力；

释放点的运动模型为：r(ξ，t)＝r_D(t).

式中，r_D(t)表示关于时间的函数；

子航天器A的运动模型为：

式中，m_A表示子航天器A的质量，F表示作用在子航天器A上的摄动力；

绞盘式释放与回收机构的运动模型为：

式中，I₁和R₁分别表示绞盘机构的转动惯量及半径，

表示系绳绞盘上所缠绕的角度，M_C表示作用在绞盘上上的控制力矩；

所述

式中，μ_e表示地球引力常数，ω表示系统的轨道角速度，R_B表示母航天器B相对地球中心的矢径；

所述释放与回收过程中满足：

所述系绳上点的应变ε和系绳的切向量τ满足：(1+ε)τ＝r′.

所述系绳中应力n满足：

式中，E表示系绳的杨氏模量，A表示系绳的截面积，α表示系绳的粘性阻尼系数；

步骤2：使用N+1个节点将系绳分为N段，第一段系绳长度为η，以后各段系绳长度为h，并对所有节点依次进行编号，计算各个节点坐标的初值：

式中，ε_i表示节点i的应变，τ_xi表示节点i处的切向量τ_i沿x轴方向的投影，τ_yi和τ_zi分别表示沿y和沿z轴的投影；

步骤3：构造系统状态微分

对于状态变量Z有：

得到函数f_z使得：

对于状态变量X有：

$\left\{\begin{array}{c}{u}_0={\left[(1+{\mathrm{\ϵ}}_0){\mathrm{\τ}}_0\right]}^{\·\·}=\frac{2}{\mathrm{\η}}g+\frac{2}{{\mathrm{\ρ\η}}^2}(n_1-n_0)+\frac{f_0+f_1}{\mathrm{\ρ\η}}-q_0\\ u_1={\left[(1+{\mathrm{\ϵ}}_1){\mathrm{\τ}}_1\right]}^{\·\·}=\frac{2}{\mathrm{\ρ\η}(h+\mathrm{\η})}(n_0-n_1)+\frac{2}{\mathrm{\ρh}(h+\mathrm{\η})}(n_2-n_1)+\frac{f_2-f_0}{\mathrm{\ρ}(h+\mathrm{\η})}-q_1\\ u_j={\left[(1+{\mathrm{\ϵ}}_j){\mathrm{\τ}}_j\right]}^{\·\·}=\frac{1}{{\mathrm{\ρh}}^2}(n_{j+1}-{2n}_j+n_{j-1})+\frac{f_{j+1}-f_{j-1}}{2\mathrm{\ρh}}-q_j(j=2,L,N-1)\\ u_N={\left[(1+{\mathrm{\ϵ}}_N){\mathrm{\τ}}_N\right]}^{\·\·}=-\frac{2}{{\mathrm{\ρh}}^2}(n_N-n_{N-1})-\frac{{2n}_N}{m_{A}h}-\frac{f_{N-1}+f_N}{\mathrm{\ρh}}+\frac{2F}{{\mathrm{hm}}_A}-q_N\end{array}\right.$

其中：

式中：

u_xj、u_yj和u_zj分别表示u_j沿x、y和z轴方向的投影，其中：

$A=\frac{1}{1+{\mathrm{\ϵ}}_j}\left[\begin{array}{cccc}{\mathrm{\τ}}_{\mathrm{yj}}^2+{\mathrm{\τ}}_{\mathrm{zj}}^2& -{\mathrm{\τ}}_{\mathrm{xj}}{\mathrm{\τ}}_{\mathrm{yj}}& -{\mathrm{\τ}}_{\mathrm{xj}}{\mathrm{\τ}}_{\mathrm{zj}}& (1+{\mathrm{\ϵ}}_j){\mathrm{\τ}}_{\mathrm{xj}}\\ -{\mathrm{\τ}}_{\mathrm{xj}}{\mathrm{\τ}}_{\mathrm{yj}}& {\mathrm{\τ}}_{\mathrm{xj}}^2+{\mathrm{\τ}}_{\mathrm{zj}}^2& -{\mathrm{\τ}}_{\mathrm{yj}}{\mathrm{\τ}}_{\mathrm{zj}}& (1+{\mathrm{\ϵ}}_j){\mathrm{\τ}}_{\mathrm{yj}}\\ -{\mathrm{\τ}}_{\mathrm{xj}}{\mathrm{\τ}}_{\mathrm{zj}}& -{\mathrm{\τ}}_{\mathrm{yj}}{\mathrm{\τ}}_{\mathrm{zj}}& {\mathrm{\τ}}_{\mathrm{xj}}^2+{\mathrm{\τ}}_{\mathrm{yj}}^2& (1+{\mathrm{\ϵ}}_j){\mathrm{\τ}}_{\mathrm{zj}}\\ (1+{\mathrm{\ϵ}}_j){\mathrm{\τ}}_{\mathrm{xj}}& (1+{\mathrm{\ϵ}}_j){\mathrm{\τ}}_{\mathrm{yj}}& (1+{\mathrm{\ϵ}}_j){\mathrm{\τ}}_{\mathrm{zj}}& -{(1+{\mathrm{\ϵ}}_j)}^2\end{array}\right]$

于是得到函数f_x和f_y使得：

将时间变量t_xz更新为t_xz+h_xz，当t_xz为h_y的整数倍时继续下一步，否则继续本步骤；

步骤4：使用四阶龙格-库塔方法求解X(t_xz+h_xz)和Z(t_xz+h_xz)：

$X(t_{\mathrm{xz}}+h_{\mathrm{xz}})=X\left(t_{\mathrm{xz}}\right)+\frac{h_{\mathrm{xz}}}{6}(k_{x1}+2k_{x2}+2k_{x3}+k_{x4})$

$Z(t_{\mathrm{xz}}+h_{\mathrm{xz}})=Z\left(t_{\mathrm{xz}}\right)+\frac{h_{\mathrm{xz}}}{6}(k_{z1}+2k_{z2}+2k_{z3}+k_{z4})$

其中：

k_x1＝h_xzf_x(t_xz，X(t_xz)，Y(t_y)，Z(t_xz))

k_z1＝h_xzf_z(t_xz，X(t_xz)，Y(t_y)，Z(t_xz))

$k_{x2}=h_{\mathrm{xz}}{f}_x(t_{\mathrm{xz}}+\frac{h_{\mathrm{xz}}}{2},X\left(t_{\mathrm{xz}}\right)+\frac{k_{x1}}{2},Y\left(t_y\right),Z\left(t_{\mathrm{xz}}\right)+\frac{k_{z1}}{2})$

$k_{z2}=h_{\mathrm{xz}}{f}_z(t_{\mathrm{xz}}+\frac{h_{\mathrm{xz}}}{2},X\left(t_{\mathrm{xz}}\right)+\frac{k_{x1}}{2},Y\left(t_y\right),Z\left(t_{\mathrm{xz}}\right)+\frac{k_{z1}}{2})$

$k_{x3}=h_{\mathrm{xz}}{f}_x(t_{\mathrm{xz}}+\frac{h_{\mathrm{xz}}}{2},X\left(t_{\mathrm{xz}}\right)+\frac{k_{x2}}{2},Y\left(t_y\right),Z\left(t_{\mathrm{xz}}\right)+\frac{k_{z2}}{2})$

$k_{z3}=h_{\mathrm{xz}}{f}_z(t_{\mathrm{xz}}+\frac{h_{\mathrm{xz}}}{2},X\left(t_{\mathrm{xz}}\right)+\frac{k_{x2}}{2},Y\left(t_y\right),Z\left(t_{\mathrm{xz}}\right)+\frac{k_{z2}}{2})$

k_x4＝h_xzf_x(t_xz+h_xz，X(t_xz)+k_x3，Y(t_y)，Z(t_xz)+k_z3)

k_z4＝h_xzf_z(t_xz+h_xz，X(t_xz)+k_x3，Y(t_y)，Z(t_xz)+k_z3).

步骤5：就将时间变量t_xz更新为t_xz+h_xz，当t_xz为h_y的整数倍时进行下一步骤，否则重复步骤4；

步骤6：采用四阶龙格-库塔方法求解Y(t_y+h_xz)：

$Y(t_y+h_y)=Y\left(t_y\right)+\frac{h_y}{6}(k_{y1}+2k_{y2}+2k_{y3}+k_{y4})$

k_y1＝h_yf_y(t，X(t_y)，Y(t_y)，Z(t_y))

$k_{y2}=h_y{f}_y(t_y+\frac{h_y}{2},X(t_y+\frac{h_y}{2}),Y\left(t_y\right)+\frac{k_{y1}}{2},Z(t_{\mathrm{xz}}+\frac{h_y}{2}))$

$k_{y3}=h_y{f}_y(t_y+\frac{h_y}{2},X(t_y+\frac{h_y}{2}),Y\left(t_y\right)+\frac{k_{y2}}{2},Z(t_y+\frac{h_y}{2}))$

k_y4＝h_yf_y(t_y+h_y，X(t_y+h_y)，Y(t_y)+k_y3，Z(t_y+h_y))

当系绳分段数N达到分段上限SupLim时，将柔性绳上编号为奇数的节点删除，使相邻分段合并；所述SupLim为偶数，范围为

当系绳分段数N达到分段下限InfLim时，则在节点i和节点i+1(i＝1，2，K，N-1)的中点处插入新的节点

其节点状态由节点节点i和节点i+1的状态取平均得到；所述InfLim范围为[1，100]；

当η≥1.5h，在节点0和节点1之间插入一个新的节点，使得新的节点与节点1之间绳段的自然长度为h，且新节点的状态由线性插值获得：

${\mathrm{\ϵ}}_{*}=\frac{h}{\mathrm{\η}}{\mathrm{\ϵ}}_0+\frac{\mathrm{\η}-h}{\mathrm{\η}}{\mathrm{\ϵ}}_1,$ ${\mathrm{\τ}}_{*}=\frac{h}{\mathrm{\η}}{\mathrm{\τ}}_0+\frac{\mathrm{\η}-h}{\mathrm{\η}}{\mathrm{\τ}}_1,$

当η＜0.5h，将第一段与第二段系绳合并起来且删去节点1；

步骤7：将时间变量t_y更新为t_y+h_y，若t_y达到仿真结束的时间则结束仿真；若未达到，则从步骤3重新开始。

Images