气囊回收系统的着陆适应性仿真分析

赵炳旗，徐振亮，吴胜宝，何欢,3，陈国平,3

（1. 机械结构力学及控制国家重点实验室，南京 210016；2. 中国运载火箭技术研究院研发中心，北京 100076；3. 南京航空航天大学 振动工程研究所；南京 210016）

半个世纪以来，人们已经成功实现了对火星、月球、小行星彗星等不同天体的探索。在探索的过程中，气囊常常被用于航天器着陆时的缓冲物。气囊作为一种弹性元件，有很好的减振吸振性能。然而，考虑到不同天体表面的复杂性，人们需要考虑较多因素，主要包括岩石尺寸、坡度大小、表面粗糙度、雷达反射率、表面承重能力等[1]。

气囊在几何上和材料上都是非线性的，在很小的外力作用下，可能引起很大的几何变形。气囊所用的薄膜材料一般为橡胶、塑胶等非线性高弹性材料[2-3]。对于气囊的计算，采用的分析方法主要有两种：一种是解析分析法，J.B.ESGAR等人从热力学方程出发，建立了气囊的解析模型，讨论了一系列气囊参数之间的关系，但并没有分析设计参数对时域缓冲特性的影响[4]。另一种是有限元仿真计算，主要的应用软件有DYTRAN和DYNA。这种有限元仿真的优势在于可对气囊的全向缓冲特性进行分析，同时可以模拟出各种着陆表面情况，缺点是不利于气囊设计参数的初始选定和设计规律的定性分析，且计算时间长[5-6]。

完全采用详细的有限元模型完成全部状态的着陆缓冲动力学仿真分析所需的计算量过于巨大。为解决该问题，文中提出了一种气囊着陆缓冲等效分析方法，将有限元仿真和理论分析相结合，借助理论分析的优点实现对气囊回收系统着陆缓冲冲击性能快速评估的目的。

1 气囊计算简化过程

如图 1所示，航天器以一定初速度撞击到斜坡上，用一个封闭的气囊对其进行减速缓冲。现将航天器等效成一个质量块m，气囊等效成一个弹簧，由于气囊刚度渐硬的特性，其弹性系数不是常量，会随着气囊压缩量的增大而增大。可利用有限元软件，建立相应的气囊模型，在无重力的条件下测出接触载荷F与气囊压缩量x的关系曲线。

斜坡表面的地形形状可通过航天器观测得到，可用一条拟合的函数曲线y=f(x)表示。航天器的着陆点具有不确定性，假设着陆点为如图2所示，以着陆点为原点，斜坡的切线方向为X轴，法线方向为轴建立局部坐标系局部坐标系与总体坐标系的夹角为θ，定义θ为撞击点处的坡度。记等效质量块在局部坐标系中的坐标为，在总体坐标系中的坐标为（x, y）。

式中：μ为斜坡表面摩擦系数。

若以气囊的中心为旋转中心，则有：

2 中心差分法

解动力学方程（1）、（2）和（3）。由于接触载荷F与气囊的压缩量有关，即F是关于的一个非线性函数，要得到方程的解析解十分困难，可以中心差分法来解这三个方程。中心差分法是一种有效的数值计算方法[7-10]，该方法的原理如下：

给定初始条件，将位移函数ut用Taylor级数展开：

解出：

现利用中心差分法计算出局部坐标系下水平和竖直方向上的位移、速度和加速度，以竖直方向为例。在轴方向上，已给出，记时间步长为在起始步中而M=m，，代入循环公式（6）中，则有

当t=2Δt时，

当t=nΔt时，

经过n+1次迭代之后，可求得局部坐标系中竖直方向质量块的位移随时间t变化的曲线，再根据式（5），可以得到速度以及加速度曲线。

3 局部坐标系与总体坐标系

由于中心差分法计算所得的位移、速度、加速度在局部坐标系中，因此需要将它们转化到总体坐标系中。如图3所示，坐标系XOY为总体坐标系，为局部坐标系。其中局部坐标系的原点在总体系下的坐标为（x0,y0） ，两个坐标系的夹角为θ，对于任意上的点，其在总体坐标系上的坐标为：

则位移转换关系式可用矩阵表示为：

同理得到速度转换关系式：

以及加速度转换关系式：

4 算例

利用 Patran软件，建立一个均压球形气囊模型，气囊半径为1 m。气囊的上方为等效质量块，质量为256 kg，气囊顶部点和上方板的中心点为同一个节点，记为点P。气囊的下方是地面，地面的节点为固定节点。气囊和质量块以初始速度20 m/s，在无重力的条件下撞击到地面上，输出节点 P的位移、速度以及加速度。以气囊的压缩量为横坐标，所受的力（质量×加速度）为纵坐标，作接触载荷-压缩量曲线关系图。根据曲线趋势，得到曲线的关系式为F=111886x2。

航天器和气囊以初始速度v0=10 m/s竖直下落，撞击到斜坡上，重力加速度为9.8 m/s2。记气囊撞击到斜坡上的时间为t0=0，建立三个方向的微分方程，编写 Matlab程序，利用中心差分法，解出微分方程（1），（2）和（3）。利用坐标转化矩阵（9），（10）和（11），得到总体坐标系下的位移、速度、加速度曲线。

输出结果以落点(0.400, 0.923)为例，此时撞击点处的斜坡坡度20.27°，如图6所示。气囊与斜坡撞击后约0.2 s离开地面，在撞击的过程中，等效质量在竖直方向的最大过载为195.4 m/s2，在撞击后，竖直方向速度的值减小，方向变为反向；等效质量在水平方向的最大过载为 66 m/s2，在撞击后，水平方向速度约为66.6 m/s。撞击过后，由于坡面的角度较大，此时转动的角速度约为0.5 rad/s。

5 仿真验证

5.1 撞击点坡度为0°

利用Patran软件，建立气囊模型，如图7所示，气囊半径为1 m。气囊的上方为等效质量块，质量为256 kg，气囊和质量块初始速度为10 m/s，方向竖直向下，重力加速度为9.8 m/s2。输出水平方向和竖直方向的过载曲线，与理论解进行对照，如图8所示。在理论解中，等效质量块水平方向的过载始终为0，竖直方向的最大过载为224.5 m/s2。在仿真解中，水平方向的过载存在一定波动，在滤波之后，过载近似为0，与理论解相当。竖直方向的过载也存在一定扰动，滤波后其最大过载约为224.5 m/s2，等于理论解。因此，当撞击点坡度为0°时，理论计算与仿真结果相一致。

5.2 撞击点坡度为20.27°

利用Patran软件，建立气囊模型，如图9所示，气囊半径为1 m。气囊的上方为等效质量块，质量为256 kg，撞击点坡度为20.27°，气囊和质量块初始速度为10号m/s，方向竖直向下，重力加速度为9.8 m/s2。输出水平方向和竖直方向的过载曲线，与理论解进行对照，如图10所示。在理论解中，等效质量块水平方向的最大过载为66.6 m/s2，竖直方向的最大过载为195.4 m/s2。在仿真解中，水平方向的过载存在一定波动，在滤波之后，最大值为 70 m/s2，与理论解的相对误差为4.8%。竖直方向的过载也存在一定扰动，滤波后其最大过载约为183 m/s2，与理论解的相对误差为 6.8%。因此，在误差允许的范围内，理论计算与仿真结果相一致。

6.3 撞击点坡度为31.24°

利用Patran软件，建立气囊模型，如图11所示，气囊半径为1 m。气囊的上方为等效质量块，质量为256 kg，撞击点坡度为31.24°，气囊和质量块初始速度为10m/s，方向竖直向下，重力加速度为9.8 m/s2。输出水平方向和竖直方向的过载曲线，与理论解进行对照，如图12所示。在理论解中，等效质量块水平方向的最大过载为92.5 m/s2，竖直方向的最大过载为157.5 m/s2。在仿真解中，水平方向的过载存在一定波动，在滤波之后，最大值为 96 m/s2，与理论解的相对误差为3.6%。竖直方向的过载也存在一定扰动，滤波后其最大过载约为148 m/s2，与理论解的相对误差为 6.4%。因此，在误差允许的范围内，理论计算与仿真结果相一致。

6 不同撞击角度情况比较

气囊和等效质量块以同一速度撞击到斜坡的不同落点时，撞击速度与坡面的夹角不同，水平方向和竖直方向的最大过载也不同。最大过载随撞击点坡度变化的曲线如图13所示。根据曲线可知，当撞击点坡度为0°时，水平方向最大过载为0，随着撞击点坡度的增大，水平方向的最大过载逐渐增大。在撞击点坡度为0°时，竖直方向最大过载的值最大，为224.5 m/s2，随着撞击点坡度增大，竖直方向的最大过载逐渐减小。此外，在不同撞击点上，比较气囊弹离地面时的角速度，如图 14所示。当撞击点的坡度在 0°到 20°之间时，气囊离地时的角速度随坡度的增加而增大；当撞击点的坡度在 20°到 30°之间时，气囊离地时的角速度随坡度的增加而增大，但增大的程度变缓。

7 结论

文中提出了一种计算气囊的理论方法。首先建立了气囊有限元模型，在无重力条件下对给定气囊施加载荷，得到接触载荷-压缩量曲线，根据曲线拟合出接触载荷与气囊压缩量的关系式。同时利用高斯函数模拟斜坡的坡面，考虑一质量块（模拟航天器）和气囊以一定初速度竖直向下撞击到该坡面上，只考虑坡度大小和表面粗糙度对气囊载荷的影响，利用中心差分法计算出质量块的位移、速度以及加速度。

将理论计算所得的结果与仿真输出的结果进行对照，在撞击点的坡度为0°，20.27°和31.24°时，得到理论的水平方向和竖直方向上的最大过载。与仿真结果对照，在误差允许的范围内，理论与仿真结果一致，从而验证了该理论计算方法的有效性。

利用该方法，分析比较不同的撞击角度下，等效质量块水平、竖直方向上最大过载以及气囊离开地面时的角速度。当撞击点坡度为 0°时，水平方向最大过载为0，随着撞击点坡度增大，水平方向的最大过载逐渐增大；竖直方向最大过载的值最大，为224.5 m/s2，随着撞击点坡度增大，竖直方向的最大过载逐渐减小；当撞击点坡度为0°时，角速度为0，气囊离开地面时的角速度逐渐增大，其增幅在 0°到20°之间较大。