0 引言

月球着陆器研发过程中需要准确地预测月面软着陆冲击载荷，从而为结构设计提供载荷输入条件。着陆器着陆冲击载荷条件的确定有两种方法。一种方法是将所有可能的输入参数的离散值进行排列组合，求出所有组合工况条件下结构所受到的最大冲击载荷作为设计载荷。采用这种方法的缺点是求出的设计载荷可能过于保守。对于着陆冲击载荷分析这样的多工况问题，这种方法存在一个难以克服的困难，即所需的仿真量过大。如果以 3 个着陆线速度、3 个着陆角速度、3 个着陆姿态角和着陆面坡度 10 个变量作为仿真输入参数，每个变量取 3 个值，可能的组合工况数高达 310个，进行这样大规模工况的仿真是不现实的。另一方法叫蒙特卡洛方法，即根据各参数的概率密度函数随机选择输入参数，根据选择的输入参数样本计算着陆冲击载荷，然后采用统计方法确定着陆冲击载荷可能的最大值，并给出相应的置信水平，这种方法可以大大降低所需的仿真工况数。蒙特卡洛法在 Surveyor，Apollo，Viking等行星探测器的设计中已经得到应用[1-7]，但由于受到当时计算速度的限制，多采用解析的刚体模型计算着陆冲击载荷。随着计算机计算速度的不断提高和有限元方法的发展成熟，使得采用有限元柔性体模型进行着陆冲击载荷的分析计算成为了可能。

 

本文通过参数化建模系统建立着陆器有限元模型，并采用具有强大非线性瞬态计算功能的 Abaqus/Explicit 作为求解器。通过与试验数据的对比验证和修正非线性有限元模型。根据

输入参数的概率密度函数随机生成 100 个着陆工况，然后利用经试验验证的有限元模型对这100 个着陆工况开展着陆冲击仿真分析，并通过 3 法则确定着陆器着陆冲击的设计载荷。

 

1 显式非线性有限元动力学分析理论

显式非线性有限元动力学模型的数值解法，是将结构在空间域进行离散，把连续的微分方程转换成有限阶的代数方程组：

 

 

 

其中 为加速度列阵， 为质量矩阵， 是单元应力引起的内部节点力列阵， 为外载荷列阵。

对求解的时间历程在时间域内离散，并采用显式积分方法进行求解。中心差分法是最常用的显式积分方法，积分步骤如下：

 

 

 

其中 为速度列阵； 为位移列阵；上标(i)表示增量步的次数，(i-1/2)和(i+1/2)是指增量步的中间值。这样，在增量步开始时提供了满足动力学平衡条件的加速度。得到了加速度后，在时间上“显式地”前推速度和位移。显式积分算法中，不需要形成整体的剪切刚度矩阵。由于是显式的前推模型的状态，所以不需要迭代和收敛准则。显式方法中可以很容易地模拟接触条件和其他一些极度不连续的情况，通过调整节点加速度来平衡在接触时的外力和内力。这些特点使得显式方法特别适合于求解着陆冲击瞬态动力学问题。虽然显式动力学求解方法的效率很高，但它却是条件稳定的，最大稳定时间积分步长 与单元的特征尺寸和波速 相关：

 

 

 

对于线弹性材料， 的计算公式如下：

 

 

 

式中： 为弹性模量； 为泊松比； 为密度。

 

2 着陆器有限元模型及参数化建模系统

2.1 着陆器有限元模型

采用 ABAQUS/CAE 建立的着陆器着陆冲击有限元模型及坐标系如图 1 所示，主要由着陆器模拟结构、四套着陆腿（包括铝蜂窝缓冲元件、拉杆缓冲元件、足垫等）和着陆面组成。着陆腿主支柱外筒及辅助支柱外筒顶端通过万向节(UJOIN)连接单元连接；辅助支柱内筒底

端与主支柱外筒之间、主支柱内筒与足垫之间通过球铰(JOIN+ROTATION)连接单元连接。着陆器足垫采用离散刚体建模。模拟结构主要由钢板和槽钢组成，钢板部分网格较粗，故采用全积分壳单元建模，而槽钢结构采用梁单元建模。着陆腿内外套筒采用壳单元建模。铝蜂窝缓冲元件采用实体单元建模，采用 Crushable Foam 体积硬化模型作为铝蜂窝缓冲器的材料模型，具体参数见文献[8]。着陆面模型由一个 7m×7m×1m 实体单元集合组成，采用Drucke-Prager/Cap 模型作为模拟月壤的本构模型，模型参数见文献。

 

 

图 1 着陆器有限元模型

Fig.1 Finite element model of lunar lander

 

2.2 着陆冲击仿真分析及试验验证

着陆仿真分析的工况见表 1。着陆器的偏航角定义为绕 X 轴转动的角度，滚转角为绕 Z轴转动的角度，俯仰角为绕 Y 轴转动的角度，三者均按右手定则方向为正。试验工况着陆器的初始角速度、俯仰角、滚转角、偏航角都近似为 0。为了与试验结果进行对比，模型中重力场取 9.8 m/s2。试验时模拟着陆器垂直方向的速度是在重力加速度作用下作自由落体运动获得的，为了节省计算时间，将模拟着陆器置于一个十分接近着陆面的位置，然后赋予初始速度来等效其在自由落体下的运动效果。着陆面模型底部和四周固支。

 

表 1 着陆工况

Tab.1 Landing scenarios

 

 

着陆缓冲装置主支柱接头载荷分析值与试验结果对比见表 2,分析值与试验实测值的最大相对误差约为 8.97%，平均相对误差约为 5.18%，分析精度满足工程要求。仿真结果与试验测量值的误差来源主要有以下几个方面：1）真实缓冲器的压缩载荷与有限元模型中设定的设计值有偏差；2）着陆腿足垫与模拟月壤之间的摩擦系数与模型中设定的固定摩擦系数存在偏差；3）力传感器的测量值本身也存在一定的测量误差；4）有限元模型做了一些理想化的假设和简化，与真实试验样机之间也存在一定的模型误差。

 

表 2 分析与试验结果对比

Tab.2 Comparison of simulation and test results

 

 

2.3 参数化建模系统

月球着陆器有限元建模涉及几何、材料和接触等非线性因素，模型比较复杂，同时由于着陆的工况很多，建模的工作量很大。为进行蒙特卡洛仿真分析，开发了月球着陆器着陆冲击参数化非线性有限元建模系统。以 Abaqus/CAE 作为建模平台，通过 Abaqus 脚本接口进行二次开发，实现探测器参数化有限元建模。

 

参数化建模系统的框架见图 2。参数化建模系统采用 Excel 2007 作为参数输入的图形用户界面（GUI），通过 VBA 宏(Macro)生成着陆器着陆工况、着陆面条件及着陆器模型的参数文件。Python 建模脚本导入建模参数文件后，通过 Python 解释器编译后进入到Abaqus/CAE 内核中执行，进入到内核中的命令将转换为 INP 文件，再经过 Abaqus/Explicit显式求解器进行分析，最后得到输出数据库（ODB）文件。Python 建模脚本是参数化建模系统的核心，它包括有数据文件、材料、零件、装配体、相互作用、连接单元、初始条件、接触、分析步、边界条件、载荷、分析作业 12 个模块。

 

 

图 2 参数化建模系统框架

Fig.2 Framework of parameterized modeling system

 

3 基于蒙特卡洛法的着陆冲击载荷分析

本节将采用蒙特卡洛方法来确定着陆器着陆冲击设计载荷。选用三个线速度(v x, v y, v z)、三个角初速度(vYaw, vPitch, vRoll)、三个姿态角(yaw, pitch, roll)和月面坡度(slope)共 10个变量作为分析的输入参数。坐标系定义见图 1。这些参数的分布概率密度函数如表 3 所示。除了偏航角和月面坡度外，其余参数均服从正态分布。偏航角(yaw)在 0~45o之间服从均匀分布，0o对应 1-2-1 着陆模式，即一个着陆腿先触地，然后 Y 轴上的两条着陆腿同时触地，+Z轴上的着陆腿最后触地；45o对应 2-2 着陆模式，即两条腿同时先触地，另两条腿最后同时着陆。由于采用模拟着陆器模型进行仿真分析，着陆器结构是关于中心轴对称的，因此偏航角在 0~45o均匀分布与在 0-360o间均匀分布是等效的。假设月面坡度服从系数为 2.113 是 3自由度的卡方(chi-square)分布，即

 

 

 

x 的概率密度函数为：

 

 

 

式中自由度 n 取 3。

 

表 3 输入参数概率密度函数

Tab.3 Probability density function of input parameters

 

 

 

根据输入参数的概率密度函数，随机生成 100 个着陆工况，提交高性能计算机进行计算。输出变量主要包括：着陆腿主支柱与着陆器结构接头处三个方向上的分力的峰值载荷 Fx，Fy，Fz；着陆腿辅助支柱与着陆器结构接头处三个分方向的分力的峰值载荷 fx，fy，fz。假设接头处三个方向的载荷峰值服从正态分布，则根据 3法则确定着陆器接头处的设计载荷。根据 3法则，设计载荷 F 可通过下式确定：

 

 

式中： 为变量的平均值； 为变量的标准差。

 

采用 3法则的前提是研究的变量服从正态分布。为了验证 100 个接头峰值载荷样本服从正态分布，通过 Kolmogorov-Smirnov 检验来验证输出的接头载荷峰值样本服从正态分布。Kolmogorov-Smirnov 检验是检验单一样本是否来自某一特定分布的方法，它是以样本数据的累积频数分布 F(x)与特定理论分布 F 0(x)比较，若两者间的差距很小，则推论该样本取自某特定分布族。即对于假设检验问题：

 

 

 

结论：当实际观测 D>D时，则接受 H1，反之则接受 H0假设。对主支柱和辅助支柱三个方向上的峰值载荷样本进行置信水平为 0.05 的 Kolmogorov- Smirnov 检验表明，6组样本均服从正态分布。主支柱接头处 X 向峰值载荷载的概率密度函数和概率分布函数的经验值和理论值见图 3 和图 4，峰值载荷样本的经验分布函数与理论正态分布函数的匹配性很好。

 

 

图 3 峰值载荷频率直方图和理论正态概率密度函数图

 

图 4 峰值载荷经验分布函数和理论正态分布函数图

 

根据 3法则得到主支柱与辅助支柱接头处峰值载荷样本均值、标准差和 3值如表 26所示, 峰值载荷超出 3值的概率仅为 0.13%。

 

表 4 接头峰值载荷均值、标准差及 3值

 

 

为确定蒙特卡洛法的有效性，对表 5 中 12 个典型着陆工况进行仿真分析。分析工况的垂直速度和水平速度均为着陆器着陆初始速度的上限。分析得到的主支柱接头 x 向峰值载荷见表 6，12 个工况中各条腿的 x 向峰值载荷最大值均小于蒙特卡洛法给出的 3值，证明了蒙特卡洛方法得到 3值是合理的。

 

表 5 典型着陆工况

Tab.5 Typical landing scenarios

 

 

表 6 着陆缓冲装置主支柱接头处 x 向峰值载荷

Tab.6 X-direction Maximum Loads at the join of main strut of langding gear

 

 

4 结束语

在通过试验数据验证分析模型正确性的基础上，利用Abaqus/CAE提供的基于Python脚本的二次开发接口，设计开发了月球着陆器着陆参数化建模系统，并利用该系数对着陆器着陆过程中的冲击载荷进行蒙特卡洛分析，通过Kolmogorov-Smirnov检验证明着陆缓冲装置接头处的峰值载荷服从正态分布，并通过3法则确定了最大着陆冲击载荷。通过对典型极限着陆工况的仿真分析证明了通过蒙特卡洛方法的有效性。

 

本文提出的基于蒙特卡洛法的着陆器着陆冲击载荷分析方法可用于确定着陆器着陆过程中使用载荷，从而为着陆器主结构的设计提供输入条件。