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

摘要

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

主权项

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+{ϵ}_0){\tau }_0\right]}^{··}=\frac{2}{\eta }g+\frac{2}{{\mathrm{\rho \eta }}^2}(n_1-n_0)+\frac{f_0+f_1}{\mathrm{\rho \eta }}-q_0\\ u_1={\left[(1+{ϵ}_1){\tau }_1\right]}^{··}=\frac{2}{\mathrm{\rho \eta }(h+\eta )}(n_0-n_1)+\frac{2}{\mathrm{\rho h}(h+\eta )}(n_2-n_1)+\frac{f_2-f_0}{\rho (h+\eta )}-q_1\\ u_j={\left[(1+{ϵ}_j){\tau }_j\right]}^{··}=\frac{1}{{\mathrm{\rho h}}^2}(n_{j+1}-{2n}_j+n_{j-1})+\frac{f_{j+1}-f_{j-1}}{2\mathrm{\rho h}}-q_j(j=2,L,N-1)\\ u_N={\left[(1+{ϵ}_N){\tau }_N\right]}^{··}=-\frac{2}{{\mathrm{\rho h}}^2}(n_N-n_{N-1})-\frac{{2n}_N}{m_{A}h}-\frac{f_{N-1}+f_N}{\mathrm{\rho 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+{ϵ}_j}\left[\begin{array}{cccc}{\tau }_{\mathrm{yj}}^2+{\tau }_{\mathrm{zj}}^2& -{\tau }_{\mathrm{xj}}{\tau }_{\mathrm{yj}}& -{\tau }_{\mathrm{xj}}{\tau }_{\mathrm{zj}}& (1+{ϵ}_j){\tau }_{\mathrm{xj}}\\ -{\tau }_{\mathrm{xj}}{\tau }_{\mathrm{yj}}& {\tau }_{\mathrm{xj}}^2+{\tau }_{\mathrm{zj}}^2& -{\tau }_{\mathrm{yj}}{\tau }_{\mathrm{zj}}& (1+{ϵ}_j){\tau }_{\mathrm{yj}}\\ -{\tau }_{\mathrm{xj}}{\tau }_{\mathrm{zj}}& -{\tau }_{\mathrm{yj}}{\tau }_{\mathrm{zj}}& {\tau }_{\mathrm{xj}}^2+{\tau }_{\mathrm{yj}}^2& (1+{ϵ}_j){\tau }_{\mathrm{zj}}\\ (1+{ϵ}_j){\tau }_{\mathrm{xj}}& (1+{ϵ}_j){\tau }_{\mathrm{yj}}& (1+{ϵ}_j){\tau }_{\mathrm{zj}}& -{(1+{ϵ}_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，且新节点的状态由线性插值获得：${ϵ}_{*}=\frac{h}{\eta }{ϵ}_0+\frac{\eta -h}{\eta }{ϵ}_1,$${\tau }_{*}=\frac{h}{\eta }{\tau }_0+\frac{\eta -h}{\eta }{\tau }_1,$当η＜0.5h，将第一段与第二段系绳合并起来且删去节点1；步骤7：将时间变量t_y更新为t_y+h_y，若t_y达到仿真结束的时间则结束仿真；若未达到，则从步骤3重新开始。