一种空间复杂柔性结构多体系统动力学建模与计算方法与流程

本发明属于空间结构动力学计算技术领域，具体涉及一种空间复杂柔性结构多体系统动力学建模与计算方法。

背景技术：

近年来，我国航天科技发展迫切需要掌握柔性空间结构在轨展开技术，以满足卫星通信、天基对地观测和深空探测等重大需求。这类空间结构一般主要由刚性杆、柔性梁、板、壳和柔性绳索等构件组成，各构件间除铰节点处的位置约束外，还存在大量同步、滑移等非线性约束。同时这类空间结构展开尺度大、柔软部件多、构形复杂，其展开过程会呈现系统大范围运动与柔软部件大变形之间的非线性耦合动力学。美国、俄罗斯等国在其航天任务中，曾发生多起结构展开失败。由于这类结构的展开动力学地面实验难度大、无法完全抵消重力影响，故其展开过程的动力学数值模拟是确保结构展开成功的关键技术。

针对此类空间复杂柔性结构多体系统，经典多体系统动力学方法的研究思路是:在柔性结构体上附着浮动坐标系,以其相对于绝对坐标系的运动来描述柔体大范围运动,通过模态缩聚近似描述柔性结构体相对浮动坐标系的变形。该类方法一般引入小转动、小变形假设,将柔性结构体的弹性变形与刚体位移分开考虑，进行解耦处理。通常这种方法只能靠大量增加单元数量来处理柔性结构体的大变形、大位移问题,从而导致计算规模大、计算效率低。此外,这种方法不能精确地反映空间复杂柔性结构展开过程中的力学本质特性。由于空间复杂柔性结构多体系统存在很多非线性约束、运动副间隙及复杂的载荷条件,目前商业软件在系统动力学建模、数值求解效率等方面存在许多困难。

技术实现要素：

有鉴于此，本发明提供了一种空间复杂柔性结构多体系统动力学建模与计算方法，为柔性结构多体系统提供了一种简单有效的动力学建模方法，提高了计算效率。

实现本发明的技术方案如下：

步骤101：建三维模型和赋予参数：根据所要求解的实物，建立该实物的三维模型，输入所建模型的材料参数、几何参数以及空间结构拓扑关系；

步骤102：网格划分：基于步骤101赋予参数后的三维模型进行网格划分，利用绝对节点坐标法、几何精确法或等几何分析方法确定网格单元类型，将划分后的三维模型定义为网格模型；

步骤103：根据步骤102划分的网格确定各个网格单元的初始节点位置和初始节点速度，基于单元初始节点位置和初始节点速度计算网格单元广义坐标向量q和广义速度向量

步骤104：根据步骤103获得的广义坐标向量和广义速度向量，设置网格模型的外部载荷及边界条件，基于网格模型的边界条件，建立网格模型的约束方程，并设置仿真总时长t以及仿真时间步长h；

步骤105：根据步骤102-104获得的网格单元类型、广义坐标向量、广义速度向量和约束方程，基于第一类拉格朗日方程建立多体系统动力学方程，如下：

其中m为系统质量阵，为广义加速度向量，f(·)为系统的弹性力向量，q(·)为系统的广义外力向量，φ(·)和φq(·)分别为步骤104获得的约束方程向量和约束方程向量对广义坐标向量的偏导数矩阵，(·)q为函数对广义坐标向量的偏导数，λ为拉格朗日乘子向量，t为仿真时间,t∈[0,t]，t为矩阵转置；

步骤106：对所建的多体系统动力学方程进行时域离散，获得时域离散后的多体系统动力学方程，对离散后的多体系统动力学方程进行迭代，并对迭代过程中生成的线性方程组进行并行求解，输出广义坐标向量q以及拉格朗日乘子向量λ，至此完成对空间复杂柔性结构多体系统动力学问题的建模与求解。

进一步地，所述步骤106具体为：

步骤601：根据步骤105得到多体系统动力学方程，对多体系统动力学方程进行离散，离散后多体系统动力学方程见公式(2)；

步骤602：对离散后的多体系统动力学方程进行求解得到第n+1步的广义坐标向量qn+1、广义速度向量和拉格朗日乘子向量λn+1；

步骤603：根据步骤602判断时间tn+1是否小于总仿真时间t，若当前时间tn+1不小于总仿真时间t，则计算停止；否则，n自动加1，重复执行步骤601-603，计算广义坐标向量qn+1、广义速度向量令时间tn+1＝tn+h，直到tn+1大于或等于t，输出广义坐标向量qn+1以及拉格朗日乘子向量λn+1，其中qn+1为q，λn+1为λ。

进一步地，对多体系统动力学方程进行离散采用如下方法：

应用广义-alpha隐式时间积分算法对多体系统动力学方程进行如下离散：

其中

式(4)中，β和γ为决定计算精度与效率的算法参数，取γ≥1/2和β≥(1/2+γ)²/4。n为迭代次数，h为仿真时间步长，为了耗散系统高频响应，引入新的算法矢量参数a，其中矢量参数a满足关系：

式(4)和式(5)中各参数的选取方法如下:

其中为算法的谱半径，αm和αf均为广义-alpha隐式时间积分算法参数。

进一步地，本步骤中用q均代表qn+1，λ均代表λn+1，所述步骤602进一步包括：

步骤6021：设置求解精度tol；

步骤6022：采用newton-rapson迭代求解所述公式(2)，获得如下线性代数方程组：

其中δq和δλ是当前时刻t系统广义坐标向量和拉格朗日乘子向量的增量，其中

和为满足如下关系的算法参数:

其中，和均为中间变量，ψ和γ均为中间矩阵，为函数对广义速度向量的偏导数；

步骤6023：根据步骤6021获得的线性方程组计算当前残差向量的模r；

其中，||·||为模；

步骤6024：根据步骤6022对线性方程组的求解更新当前广义坐标向量以及拉格朗日乘子向量；

步骤6025：根据步骤6023获得的残差向量模判断迭代结果是否收敛，若r>tol，则根据步骤6024更新获得的广义坐标向量以及拉格朗日乘子向量，重复执行步骤6022-6025，直至r<tol或r＝tol，输出步骤6024获得的广义坐标向量q以及拉格朗日乘子向量λ。

有益效果：

1)本发明采用绝对节点坐标法、几何精确梁理论和等几何分析等非线性有限元法对空间复杂柔性结构进行离散，充分考虑了柔性结构的运动与变形耦合问题，为空间复杂柔性结构多体系统的动力学建模提供了手段。

2)本发明采用区域分解和并行计算技术，大幅度提高了空间复杂柔性结构多体系统动力学仿真计算效率。

附图说明

图1为对含复杂多柔体系统动力学问题进行建模与计算方法流程图。

图2为某桁架索网天线结构示意图。

图3为本发明所开发的部分可选的单元类型示意图。

图4为轴承轴瓦模型结构示意图。

图5为同步齿轮模型结构示意图。

图6为移动副模型结构示意图。

图7为对系统动力学方程高效数值积分求解算法流程图。

图8为对系统离散非线性方程高效求解算法流程图。

图9为本发明建模和求解算法方法应用于某桁架索网天线，并进行数值仿真天线展开部分时刻结构示意图。

具体实施方式

下面结合附图并举实施例，对本发明进行详细描述。

对于空间复杂柔性结构多体系统，通常由众多柔性杆件和柔性索网组成，在展开过程中这些构件会产生很大变形和大范围转动，此时基于小变形和小转动假设的多体建模分析方法已不能满足工程需求。对于此类柔性构件的建模主要采用非线性有限元法：如绝对节点坐标法、几何精确法和等几何分析方法等。其中shabana等基于连续介质力学和有限元法提出了绝对节点坐标法(ancf)。与经典多体系统动力学方法相比，该方法在惯性坐标系下采用统一的插值函数描述柔体的刚体运动与弹性变形，非常适用于描述同时具有大转动和大变形的柔体动力学。此外，几何精确梁理论是simo等以有限应变梁理论为基础发展的一种非线性有限元方法，该法能够高效并准确地处理大变形与大转动问题，近年来也被广泛应用于空间可展开结构的动力学仿真。为了统一几何模型与力学分析模型，实现两者的无缝集成，在有限元领域，hughes等提出的一种新的以样条理论为基础的数值计算方法：等几何分析(iga)。

综上所述，在空间复杂柔性结构多体系统动力学展开分析阶段，如何考虑并解决柔性构件变形与大范围运动的耦合问题，并提高数值仿真计算效率成为航天科技发展亟待解决的技术问题。本发明提出了一种对含大变形构件的空间结构多体系统动力学建模与计算方法。

如图1所示，实施典型案例为对如图2所示的某环形桁架索网天线机构展开动力学建模与计算方法。

该方法包括如下步骤：

步骤101：建三维模型和赋予参数：根据所要求解的实物，建立该实物的三维模型，输入所建模型的材料参数、几何参数以及空间结构拓扑关系。具体包括：根据环形桁架索网天线的几何参数以及空间结构拓扑关系，在proe软件中建立环形桁架索网天线三维模型，并输入环形桁架天线不同部件的材料参数。由于金属丝反射网密度很小，且附着在前张力网上，因而忽略金属丝反射网的影响。

步骤102：网格划分：基于步骤101赋予参数后的三维模型进行网格划分，利用绝对节点坐标法、几何精确法或等几何分析方法确定网格单元类型，将划分后的三维模型定义为网格模型；根据模型性质以及精度要求，选取合适的单元类型，如图3所示，刚体模型有：基于自然坐标方法描述的刚体单元、基于自然坐标方法描述的缩减刚体单元，基于绝对节点坐标方法描述的参考节点刚体单元；柔体单元有：基于绝对节点坐标方法描述的全参数梁单元、基于绝对节点坐标方法描述的缩减梁单元、基于绝对节点坐标方法描述的斜率不连续的梁单元、基于绝对节点坐标描述的矩形薄板单元以及基于绝对节点坐标描述的三角形薄板单元等，其中绝对节点坐标法、几何精确法或等几何分析方法均具有解决耦合特性的特点。

步骤103：确定模型广义坐标向量和广义速度向量：根据步骤102划分的网格确定各个网格单元的初始节点位置和初始节点速度，基于单元初始节点位置和初始节点速度计算网格单元广义坐标向量q和广义速度向量通过将环形桁架索网天线收拢到一定角度，从而确定索网和桁架初始广义坐标和广义速度。

步骤104：根据步骤103获得的广义坐标向量和广义速度向量，设置网格模型的外部载荷及边界条件，基于网格模型的边界条件，建立网格模型的约束方程，并设置仿真总时长t以及仿真时间步长h。该环形桁架索网天线仿真总时长为500s，仿真时间步长设为1e-4s。其中外部载荷为索网的预张力，约束为索网和桁架之间的铰接、桁架杆件之间的运动副以及桁架与伸展臂之间的固接。下面简要介绍三种主要运动副的建模：旋转铰约束、同步齿轮约束和移动副约束。

其中关节旋转铰约束建立方式如下：对旋转的关节处建立刚性的轴承轴瓦模型，如图4所示。刚性的轴承与轴瓦与相连接的构件之间为固定连接，则相连构件的旋转可以通过轴承轴瓦的相对转动来描述。采用自然坐标描述的轴承轴瓦之间的旋转铰约束方程可表示为：

其中φ¹为旋转铰约束，分别由位移约束φ1和旋转角约束φ2组成。和分别为全局坐标系下轴承和轴瓦端点位置坐标。e¹和e²分别是局部坐标系下轴瓦和轴承的单位向量，θ为轴承绕轴瓦的转动角度，其中

其中lb和lj分别是轴承和轴瓦的长度，u^b和v^b是位于轴瓦端点b的单位向量，u^j和v^j是位于轴承端点j的单位向量。

同步齿轮约束建立方式如下：

如图5所示，同步齿轮约束保证相邻两个面的展开角度一致，即第一跨的横杆(ab杆)与竖杆之间的夹角和第二跨的横杆(cd杆)与竖杆之间的夹角相等，其具体的约束方程见下式:

其中φ²为同步齿轮约束，和分别是竖杆与两横杆的轴线方向矢量。

如图6所示，移动副约束建立如下：移动副约束需保证细杆gh能在粗杆ef中自由滑动，同时细杆与粗杆在桁架展开任意位置时均在四边形abcd内。其具体的约束方程见下式：

其中φ³为移动副约束。φ3保证e点与a点重合，c点与g点重合，其中r^a、r^e、r^c和r^g分别为点a、e、c和g在全局坐标系下的位置坐标。φ4保证两斜杆能相互滑动，其中和分别为点e、h、f和g沿斜杆方向矢量，l为斜杆长度。φ4中的最后一项是为了保证两斜杆在平面abcd内，这样式φ3与φ4就构成了移动副的约束方程，至此，φ¹、φ²和φ³构成网络模型的约束方程向量。

步骤105：根据步骤102-104获得的网格单元类型、广义坐标向量、广义速度向量和约束方程，基于第一类拉格朗日方程建立多体系统动力学方程，如下：

其中m为系统质量阵，为广义加速度向量，是广义坐标向量对时间的二阶导数，f(·)为系统的弹性力向量，q(·)为系统的广义外力向量，φ(·)和φq(·)分别为步骤104获得的网络模型的约束方程向量和约束方程向量对广义坐标向量的偏导数矩阵，(·)q为函数对广义坐标向量的偏导数，λ为拉格朗日乘子向量，t为仿真时间,t∈[0,t]，t为矩阵转置。

步骤106：利用广义-alpha隐式时间积分算法，对所建的多体系统动力学方程进行时域离散，获得时域离散后的多体系统动力学方程，采用newton-rapson迭代对离散后的多体系统动力学方程进行迭代，并对迭代过程中生成的线性方程组进行并行求解，输出广义坐标向量q以及拉格朗日乘子向量λ。

如图7所示，步骤601-603描述了所述步骤106，具体为：

步骤601：根据步骤105得到多体系统动力学方程，应用广义-alpha隐式时间积分算法进行如下离散：

其中

式(7)中，β和γ为决定计算精度与效率的算法参数，取γ≥1/2和β≥(1/2+γ)²/4。n为迭代次数，h为时间积分步长，为了耗散系统高频响应，引入新的算法矢量参数a，其中矢量参数a满足关系：

式(7)和式(8)中各参数的选取方法如下:

其中为算法的谱半径，决定着算法能量耗散分布的频率范围，通常取αm和αf均为广义-alpha隐式时间积分算法参数。

步骤602：将式(7)带入式(6)得到离散形式的非线性方程组，并对非线性方程组进行求解得到广义坐标向量qn+1、广义速度向量和拉格朗日乘子向量λ；

步骤603：根据步骤602判断时间tn+1是否小于总仿真时间t，若当前时间tn+1不小于总仿真时间t，则计算停止；否则，n自动加1，令时间tn+1＝tn+h重复执行步骤601-603，计算广义坐标向量qn+1、广义速度向量直到tn+1大于或等于t。输出广义坐标向量qn+1以及拉格朗日乘子向量λn+1，其中qn+1为q，λn+1为λ。

如图8所示，步骤6021-6025进一步地描述了步骤602，具体为：

步骤6021：设置求解精度tol；

步骤6022：采用newton-rapson迭代求解公式(6)，获得如下线性代数方程组：

其中δq和δλ是当前时刻t系统广义坐标向量和拉格朗日乘子向量的增量，其中

和为满足如下关系的算法参数:

其中，和均为中间变量，ψ和γ均为中间矩阵，为函数对广义速度向量的偏导数；

步骤6023：根据步骤6021获得的线性方程组计算当前残差向量的模r；

其中，||·||为模；

步骤6024：根据步骤6022对线性方程组的求解更新当前广义坐标向量以及拉格朗日乘子向量；

步骤6025：根据步骤6023获得的残差向量模判断迭代结果是否收敛，若r>tol，则根据步骤6024更新获得的广义坐标向量以及拉格朗日乘子向量，重复执行步骤6022-6025，直至r<tol或r＝tol，输出步骤6024获得的广义坐标向量q以及拉格朗日乘子向量λ，并存储。至此完成对某环形桁架索网天线机构展开动力学建模与高效计算。图9给出了某环形桁架天线展开过程中6个指定时刻的系统构型，该图明显展示了环形桁架天线展开过程中的不同步现象相应的地面试验结果验证了该数值仿真结果的有效性。

综上所述，以上仅为本发明的较佳实施例而已，并非用于限定本发明的保护范围。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。