文献[3-4]分别针对单组和多组着陆工况，以提升着陆器的软着陆性能为目标优化了着陆器的缓冲机构。文献[5]设计了着陆器用液压式缓冲器，并提出了缓冲器的自适应控制方法，且通过试验装置验证了缓冲器的性能。文献[6-8]基于磁流变阻尼器设计了着陆器的变阻尼缓冲器及其控制策略，并通过理论分析与试验验证的方法总结了磁流变液变阻尼缓冲器的优点。文献[9-10]将磁流变液变阻尼缓冲器引入着陆器整机之中，通过仿真计算的方法验证了变阻尼缓冲器可提升着陆器的软着陆性能。以上文献的研究重点多放在变阻尼缓冲器的设计、控制方法与性能验证方面，对变阻尼缓冲器的缓冲特性参数进行优化的研究鲜有文献报道，而为了获得性能最佳的变阻尼缓冲器，有必要对其缓冲特性参数进行优化。
有鉴于此，本文以参数可控的变阻尼缓冲器为研究对象，将其应用于着陆器着陆缓冲机构的设计中。通过对着陆器在不确定着陆工况下软着陆性能的仿真分析，验证了缓冲器的缓冲性能满足要求。而后，为提高缓冲器的应用性能，结合响应面法与蒙特卡罗法，研究了变阻尼缓冲器的缓冲特性参数优化方法，并加以实施。
1 变阻尼缓冲器及其控制策略 目前有多种方案可实现阻尼调节，考虑技术成熟度，本文基于典型磁流变阻尼器设计了一种腿式着陆器用缓冲器，其结构如图 1所示[8-10]。这种缓冲器相比于机械式缓冲器的优点为可通过半自动控制技术实现自适应控制，以提高着陆器的地形适应能力。变阻尼缓冲器的主要缓冲元件为弹簧和磁流变液，其中弹簧提供恢复力，磁流变液提供阻尼力，通过控制电流强度来控制线圈的磁场强度，从而控制阻尼器的阻尼系数，实现缓冲力大小的调节。
 |
| 图 1 变阻尼缓冲器示意图[8-10] Fig. 1 Schematic diagram of variable damping buffer[8-10] |
| 图选项 |
设定变阻尼缓冲器作为装有4套缓冲机构的着陆器的主缓冲器，为了实现变阻尼缓冲器的控制，定义如图 2所示的坐标系[11]。
 |
| 图 2 坐标系定义 Fig. 2 Definition of coordinate systems |
| 图选项 |
图 2中，OS-XSYSZS为与着陆面固联的全局坐标系；OC-XCYCZC为质心控制坐标系。同时，设全局坐标系依次绕ZS-XS-YS轴旋转到与质心控制坐标系平行位置时绕各坐标轴的旋转角分别为φ(绕XS轴的转角)、θ(绕YS轴的转角)和ψ(绕ZS轴的转角)。
根据文献[10-11]的研究成果，变阻尼缓冲器的缓冲力模型可表示为
 | (1) |
式中：fs为缓冲器缓冲力；c为缓冲器等效阻尼系数；k为弹簧刚度系数；d为缓冲器的相对缓冲行程。
为了确保着陆器的软着陆缓冲性能及控制的灵活性和及时性，本文采用基于能量最低原则的跳跃控制策略实现缓冲器阻尼系数的控制[10]。针对常见的缓冲机构对称布置的着陆器，阻尼系数控制函数如式(2)所示。
 | (2) |
式中：ci为第i个缓冲器的阻尼系数；cmin和cmax分别为可变阻尼系数的最小值和最大值；di为第i个缓冲器的相对缓冲行程。
2 着陆器软着陆性能判据及工况 2.1 软着陆性能判据 为了研究着陆器的软着陆性能，首先需要确立软着陆性能的判据。根据常见的深空探测任务要求，着陆器主要软着陆性能判据为以下4点：
1) 着陆器具有较强的抗翻倒能力。即着陆过程中着陆器质心距含任意相邻两足垫中心点的竖直平面的最小距离(T)应大于0，否则认为着陆器翻倒，且T越大，着陆器抗翻倒能力越强[12]。
2) 加速度过载特性。考虑着陆器搭载的探测设备承受过载的能力有限，为保证探测任务的顺利进行，软着陆过程中的最大加速度过载(L)不应超过许用值，且L值越小越好。
3) 着陆器底部具有较好的抗损坏能力。即着陆过程中，位于着陆器主体结构底部的尾喷管的底面中心点距着陆面的最小距离(U)应大于许用值。且U越大，着陆器底面越不易受损。
4) 若采用变阻尼缓冲器作为整个着陆器的主要吸能装置，需要保证其最大缓冲行程(S)小于许用值，否则缓冲器可能发生损坏。
2.2 软着陆工况的确定 设着陆器采用非零高度关机着陆模式实现软着陆，其首先缓速运动到距离着陆面一定高度处，若地貌识别系统判断该地点满足着陆条件，则发动机关闭，着陆器自由降落至着陆面。不考虑着陆器与着陆面接触之前的主动控制误差，着陆器与着陆面接触时水平速度、三轴角速度为0、竖直速度为vx。考虑着陆器装有地貌识别系统，足垫落入凹坑或撞击石块的可能性极小。因此，可将着陆工况参数简化为着陆面摩擦系数(μ)、着陆面坡度(α)、着陆偏航角(θp)如图 3所示。
 |
| 图 3 着陆工况示意图 Fig. 3 Schematic diagram of landing conditions |
| 图选项 |
考虑着陆器中心对称，θp取值在0~45°范围内即包含评判软着陆性能所需的全部工况。
3 着陆器软着陆性能分析 3.1 着陆器整机动力学建模 本文以典型的悬臂腿式着陆器为研究对象，验证变阻尼缓冲器的缓冲性能，着陆器整机构型如图 4所示，其主要由主体结构和缓冲机构组成，缓冲机构共4套，构型和尺寸均相同，均匀分布在主体结构周围。缓冲机构由主支柱、辅助支柱和足垫组成，主支柱的上端通过万向节与主体固连，主支柱的下端通过球铰与足垫连接。辅助支柱通过万向节与主体相连接，下端通过球铰与主支柱相连接。其中，主支柱为变阻尼缓冲器，辅助支柱内填充铝蜂窝以进行辅助缓冲。
 |
| 图 4 着陆器整机构型 Fig. 4 Complete machine configuration of lander |
| 图选项 |
本文应用ADAMS软件实现着陆器动力学仿真模型的建立，仿真模型中的重力加速度设置为月球重力加速度。为了提高仿真效率并方便与2.2节中的工况参数结合，动力学仿真模型从着陆器首个足垫与着陆面接触的瞬间开始进行仿真计算。此外，由于主要关注各构件间相对运动状态及整体姿态变化，建模时将着陆器的各构件均简化为刚体。足垫与着陆面之间的作用力采用接触力模型模拟，其相互作用力可分解为法向碰撞力与切向摩擦力，法向碰撞力利用非线性弹簧阻尼模型模拟，切向摩擦力采用库仑摩擦模型模拟[4, 13]。
针对第1节的控制策略，在ADAMS中实时采集着陆器着陆过程中当前YS和ZS轴方向的角速度
3.2 软着陆动力学分析 参考文献[10]的相关研究成果，确定变阻尼缓冲器缓冲特性参数的取值如表 1所示。此外，本文算例考虑的着陆工况参数取值及其概率分布规律[14-16]如表 2所示。同时，设定着陆器整机质量为1 710 kg，vx取3.4 m/s，U的许用值≥350 mm，S的许用值≤150 mm，L的许用值≤13g(g=9.8 m/s2)。
表 1 变阻尼缓冲器参数 Table 1 Parameters of variable damping buffer
| 参数 | 取值 | 取值范围 |
| k/(kN·m-1) | 50 | 30~ 70 |
| cmax/(kN·s·m-1) | 50 | 30~ 70 |
| r | 0.5 | 0.3~0.7 |
| 注：r=cmin/cmax—阻尼系数比。 | ||
表选项
表 2 着陆工况参数取值 Table 2 Value of of landing condition parameters
| 参数 | 分布规律 | 分布参数 | (μ0-6σ，μ0+6σ] |
| α/(°) | 正态分布 | μ0=6，σ=0.667 | (2，10] |
| μ | 正态分布 | μ0=0.4，σ=0.016 7 | (0.3，0.5] |
| θp/(°) | 均匀分布 | [0, 45] | |
| 注：μ0—均值；σ—标准差。 | |||
表选项
基于上述ADAMS模型，本文采用蒙特卡罗法对着陆器在不确定条件下的着陆工况的软着陆性能进行分析[17]。具体分析流程为：首先根据表 2中的参数分布采用描述性采样方法抽取1 000个样本点，然后代入上述仿真模型进行计算，得到软着陆性能指标值的均值、标准差与超出许用值的概率分布。通过仿真计算，得到软着陆性能分析结果如表 3所示。
表 3 软着陆性能分析结果 Table 3 Analysis results of soft landing performace
| 性能指标 | 均值 | 标准差 |
| U/mm | 436.281 | 1.207 |
| S/mm | 53.555 | 2.801 |
| L/g | 9.625 | 0.290 |
| T/mm | 1 703.846 | 37.111 |
| 注：超出许用值概率均为0。 | ||
表选项
仿真结果显示，将变阻尼缓冲器引入着陆器中具有可行性。由于变阻尼缓冲器初始缓冲特性参数的选取具有较大的经验性，为了获得软着陆性能最佳的着陆器，下面将对变阻尼缓冲特性器的缓冲参数进行优化。
4 变阻尼缓冲器缓冲特性参数优化 基于3.2节的分析，本文选取k、cmax和r作为变阻尼缓冲器的设计变量对其进行参数优化，以提升着陆器的软着陆性能。
4.1 响应面模型的建立 本文的优化过程中，将继续考虑着陆工况的不确定性。此时，若仍利用ADAMS模型结合蒙特卡罗法进行迭代分析将消耗大量的时间成本。故本文建立响应面模型来模拟U、S、L和T4个输出量与k、cmax、r、α、μ和θp 6个输入变量的映射关系，结合响应面模型与蒙特卡罗法将有效提高研究着陆器软着陆性能的效率[11]。本文应用包含6个输入变量的不完全三阶多项式作为近似函数，其表达式为
 | (3) |
式中：x为输入变量向量；xi为输入变量；n为输入变量的个数；βi为多项式的各项系数。抽取n个输入变量的样本点，并将样本点代入动力学仿真模型计算出对应的输出值。将样本点的输入、输出值代入式(3)，得到响应面模型的向量表达式为
 | (4) |
式中：Y为利用动力学仿真模型计算出的输出向量；X为输入变量矩阵；?为近似函数值向量；ε为近似函数误差向量；β为近似函数的系数向量。利用最小二乘法拟合可得到β的数值，从而获得具体的近似函数表达式[18-19]。
根据6σ原则，可确定表 3所示的各着陆工况参数的取值范围，利用优化拉丁超立方实验设计在k、cmax、r、α、μ和θp 6个输入量的取值范围内抽取35组输入变量样本点，代入ADASM模型，对应计算出U、S、L和T。基于仿真结果，可拟合出U、S、L和T的响应面近似函数。本文在拟合响应面近似函数的过程中，进行了关键项筛选，以便去掉对响应面近似函数的精度产生负面影响的项，提升响应面模型的拟合精度。
本文采用均方根相对误差值(RMSE)和决定系数(R2)来检验响应面模型对样本数据的模拟精度，分别为
 | (5) |
 | (6) |
式中：nS为样本点总数；yi为输入样本点计算出的输出值；?i为响应面模型估计值；y为yi的平均值。RMSE趋于0，R2趋于1说明响应面精度较高。响应面拟合精度检验结果如表 4所示，结果显示本文所用响应面模型可较准确地描述k、cmax、γ、α、μ和θp与U、S、L和T之间的映射关系。
表 4 响应面模型精度分析 Table 4 Accuracy analysis of response surface model
| U | 1.525 | 0.998 |
| S | 10.25 | 0.999 |
| L | 1.856 | 0.999 |
| T | 1.321 | 0.999 |
表选项
利用响应面模型结合蒙特卡罗法分析着陆器软着陆性能的过程为：将k、cmax和r设置为定值，此时响应面模型的变量仅为α、μ和θp；再按表 2中α、μ和θp的分布规律抽取一定数量的样本点，代入响应面模型计算；最后根据计算结果可得到U、S、L和T的均值、标准差等统计学特征。
4.2 缓冲特性参数优化设计 为了提升着陆器在不确定工况下的软着陆性能，本文以最大化U、T的均值μU、μT，最小化L的均值μL与U、T和L的标准差σU、σT和σL为优化目标，以蒙特卡罗模拟过程中U、T的最小取值与S、L的最大取值不超过许用值作为约束条件建立优化数学模型[20]。其中，性能指标值的均值与标准差通过加权处理整合为一个目标函数，加权系数ζ可以根据具体设计要求做出调整。此外，为了防止优化后着陆器的性能指标值与许用边界重合，优化过程中，将U的许用值提升为360 mm，T的许用值提升为100 mm，S的许用值降为140 mm，L的许用值降为12 g。综上，本文取加权系数ζ=0.2，建立优化数学模型如下：
 | (7) |
式中：xO=(k，cmax，r)T为设计变量，xO(L)与xO(U)分别表示设计变量的下界与上界；Umin、Tmin与Smax、Lmax分别为蒙特卡罗模拟过程中U、T的最小值与S、L的最大值。本文应用第二代非劣排序遗传算法(NSGA-Ⅱ)进行迭代计算，算法参数配置如表 5所示。
表 5 优化参数设置 Table 5 Value of optimal parameters
| 参数 | 种群数 | 进化代数 | 交叉指数 | 变异指数 | 交叉概率 |
| 数值 | 12 | 20 | 10 | 20 | 0.9 |
表选项
每轮迭代计算均利用3.1节提到的响应面法进行100 000次蒙特卡罗模拟，输出U、L和T的均值与标准差，同时输出U、T的最小值与S、L的最大值。经过迭代计算，得到该优化问题的Pareto最优解集，并绘制Pareto前沿曲面如图 5所示。
 |
| 图 5 Pareto前沿曲面 Fig. 5 Surface of Pareto front |
| 图选项 |
在Pareto最优解集中，选取使着陆器加速度过载最小的解作为最优解，即xO*=(k*，cmax*，r*)T=(57 460.663 N/m，42 509.866 N·s/m，0.638)T。将该组解代入仿真模型，再次抽取1 000个样本点进行蒙特卡罗模拟。通过仿真计算，得到软着陆性能分析结果如表 6所示。
表 6 优化后软着陆性能分析结果 Table 6 Analysis results of soft landing performace after optimizaiton
| 性能指标 | 均值 | 标准差 |
| U/mm | 437.505 | 1.284 |
| S/mm | 49.588 | 2.564 |
| L/g | 9.485 | 0.234 |
| T/mm | 1 701.452 | 38.067 |
| 注：超出许用值概率均为0。 | ||
表选项
根据优化前后的2 000组蒙特卡罗模拟结果，做出U、L、T的散点图如图 6所示。
 |
| 图 6 U-L、U-T及T-L散点图 Fig. 6 Scatter plot of U-L, U-T and T-L |
| 图选项 |
对比表 3与表 6的结果可知：优化后，在U、S、L和T超出许用值概率依旧为0的前提下，U的均值提升了0.281%、标准差增加了6.379%，L的均值降低了1.455%、标准差降低了19.310%，T的均值降低了0.141%、标准差增加了2.576%。
由图 6可直观看出：优化后，U的标准差虽然有所增加，但其取值更加远离许用边界，着陆器底面抗损坏能力增强；优化后，L取值更加远离许用边界，着陆器的加速度过载得到削弱；优化后，T的均值虽然有所减小且标准差有所增加，但T与许用边界的距离基本保持不变；从而可知，优化后着陆器的软着陆综合性能得到提升。
5 结论 1) 本文以一种参数可控的腿式着陆器用变阻尼缓冲器为研究对象，应用ADAMS软件建立了装有变阻尼缓冲器的着陆器动力学仿真模型，并结合仿真模型与蒙特卡罗法验证了变阻尼缓冲器应用于着陆器中的可行性。
2) 基于仿真模型计算结果，构造了描述变阻尼缓冲器的缓冲特性参数、着陆工况参数和软着陆性能指标值之间复杂关系的响应面模型，并利用响应面模型、蒙特卡罗法与NSGA-Ⅱ优化了变阻尼缓冲器的相关参数。利用仿真模型验证，优化后的变阻尼缓冲器使着陆器的软着陆综合性能得到提升。
3) 本文的研究内容对着陆器变阻尼缓冲器的设计具有借鉴意义，且本文的研究方法具有通用性，可用于其他机械系统的优化设计中。
参考文献
| [1] | 李萌.腿式着陆缓冲装置吸能特性及软着陆过程动力学仿真研究[D].哈尔滨: 哈尔滨工业大学, 2013: 1-10. LI M.Research on energy absorbers of legged-type lander and dynamic simulation on its soft landing process[D].Harbin: Harbin Institute of Technology, 2013: 1-10(in Chinese). http://cdmd.cnki.com.cn/Article/CDMD-10213-1013045385.htm |
| [2] | 杨建中, 曾福明, 满剑锋, 等. 嫦娥三号着陆器着陆缓冲系统设计与验证[J]. 中国科学:技术科学, 2014, 44(5): 440-449. YANG J Z, ZENG F M, MAN J F, et al. Design and verification of the landing impact attenuation system for Chang'E-3 lander[J]. Scientia Sinica Technologica, 2014, 44(5): 440-449. (in Chinese) |
| [3] | 王家俊, 王春洁, 宋顺广. 基于响应面法的月球着陆器软着陆性能优化[J]. 北京航空航天大学学报, 2014, 40(5): 707-711. WANG J J, WANG C J, SONG S G. Performance optimization of lunar lander based on response surface methodology[J]. Journal of Beijing University of Aeronautics and Astronautics, 2014, 40(5): 707-711. (in Chinese) |
| [4] | 吴宏宇, 王春洁, 丁建中, 等. 基于多工况的新型着陆器软着陆性能优化[J]. 北京航空航天大学学报, 2017, 43(4): 776-781. WU H Y, WANG C J, DING J Z, et al. Soft landing performance optimization for novel lander based on multiple working conditions[J]. Journal of Beijing University of Aeronautics and Astronautics, 2017, 43(4): 776-781. (in Chinese) |
| [5] | LUCAS G H, ROBERT H D, VELORIA J M. Modeling and validation of a navy A6-Intruder actively controlled landing gear system:NASA/TP-1999-209124[M]. Washington, D.C.: NASA, 1999. |
| [6] | MIKULOWSKI G M, HOLNICKISZULC J. Adaptive landing gear concept-feedback control validation[J]. Smart Materials and Structures, 2007, 16(6): 2146-2158. DOI:10.1088/0964-1726/16/6/017 |
| [7] | MIKULOWSKI G, JANKOWSKI L. Adaptive landing gear:Optimum control strategy and potential for improvement[J]. Shock and Vibration, 2015, 16(2): 175-194. |
| [8] | 牛伯瑶.月球探测器着陆系统磁流变缓冲器的结构设计与优化[D].重庆: 重庆大学, 2016: 23-55. NIU B Y.The structure design and optimization of magnetorheological buffer of lunar landing system[D].Chongqing: Chongqing University, 2016: 23-55(in Chinese). http://www.wanfangdata.com.cn/details/detail.do?_type=degree&id=D01014719 |
| [9] | 汪岸柳.月球着陆器软着陆动力学与半主动控制研究[D].南京: 南京航空航天大学, 2011: 80-109. WANG A L.Investigation on the soft-landing dynamics and semi-active control for lunar lander[D].Nanjing: Nanjing University of Aeronautics and Astronautics, 2011: 80-109(in Chinese). http://www.wanfangdata.com.cn/details/detail.do?_type=degree&id=Y2260950 |
| [10] | MAEDA T, OTSUKI M, HASHIMOTO T, et al. Attitude stabilization for lunar and planetary lander with variable damper[J]. Journal of Guidance, Control, and Dynamics, 2016, 39(8): 1790-1804. DOI:10.2514/1.G000325 |
| [11] | 吴宏宇, 王春洁, 丁宗茂, 等. 两种着陆模式下的着陆器缓冲机构构型优化[J]. 宇航学报, 2017, 38(10): 1032-1040. WU H Y, WANG C J, DING Z M, et al. Configuration optimization of landing gear under two kinds of landing modes[J]. Journal of Astronautics, 2017, 38(10): 1032-1040. (in Chinese) |
| [12] | ZUPP G A, DOIRON H H.A mathematical procedure for predicting the touchdown dynamics of a soft landing vehicle: NASA-TN-D-7045[R].Washington, D.C.: NASA, 1971. |
| [13] | NOHMI M, MIYAHARA A.Modeling for lunar lander by mechanical dynamics software: AIAA-2005-6416[R].Reston: AIAA, 2005. |
| [14] | MERCHANT D H, SAWDY D T. Monte Carlo dynamic analysis for lunar module landing loads[J]. Journal of Spacecraft and Rockets, 1971, 8(1): 48-55. DOI:10.2514/3.30216 |
| [15] | MURACA R J, CAMPBELL J W, KING C A.A Monte Carlo analysis of the Viking lander dynamics at touchdown: NASA-TN-D-7959[R].Washington, D.C.: NASA, 1975. |
| [16] | 宋顺广, 王春洁. 基于蒙特卡罗法的月球探测器着陆稳定性分析[J]. 北京航空航天大学学报, 2013, 39(9): 1348-1352. SONG S G, WANG C J. Landing stability analysis of the lunar lander based on Monte Carlo approach[J]. Journal of Beijing University of Aeronautics and Astronautics, 2013, 39(9): 1348-1352. (in Chinese) |
| [17] | 李铁柱, 李光耀, 陈涛, 等. 基于Kriging近似模型的汽车乘员约束系统稳健性设计[J]. 机械工程学报, 2010, 46(22): 123-129. LI T Z, LI G Y, CHEN T, et al. Robustness design of occupant restraint system based on Kriging model[J]. Journal of Mechanical Engineering, 2010, 46(22): 123-129. (in Chinese) |
| [18] | KURTARAN H, ESKANDARIAN A, MARZOUGUI D. Crashworthiness design optimization using successive response surface approximations[J]. Computational Mechanics, 2002, 29(4-5): 409-421. DOI:10.1007/s00466-002-0351-x |
| [19] | LI M, DENG Z Q, GUO H W, et al. Optimizing crashworthiness design of square honeycomb structure[J]. Journal of Central South University of Technology, 2014, 21(3): 912-919. DOI:10.1007/s11771-014-2018-0 |
| [20] | TSUIK L. An overview of Taguchi method and newly developed statistical methods for robust design[J]. ⅡE Transactions, 1992, 24(5): 44-57. |
