飞机起落架的动力学分析与拓扑优化研究

赵知辛，王 琨，汪 杰，张昌明

（1.陕西理工大学机械工程学院，陕西 汉中723000；2.陕西省工业自动化重点实验室，陕西 汉中723000）

1 引言

飞机起落架是实现飞机着陆、滑跑以及降落的重要装置，并在起降、滑跑过程中用于支撑飞机的重量［1，2］。因此对起落架进行落震仿真分析与结构优化设计尤为重要。在飞机飞行过程中，起落架是收在飞机机身内部的，如果起落架的体积很大，意味着收放空间就要越大，这就会造成起落架的收放困难［3，4］。并且随着起落架的重量增加，对燃油量的消耗也会增多，实验数据表明，每减少10kg的飞机重量，就能够节省3925kg的燃油量并少排放4吨CO2。因此对起落架进行减重是当前的研究热点。国内外已经有学者对飞机起落架进行研究，在国外，文献［5］提出一种轻型紧凑型直升机着陆支柱的优化设计方法。文献［6］提出了一种考虑强度、稳定性和挠度限制最小重量的飞机壁板层复合材料结构设计方法。文献［7］提出了一种解决航天设计问题的拓扑优化算法。在国内，文献［8］建立了飞机起落架带锁撑杆支杆的三维模型，利用Matlab对支杆的截面尺寸进行了优化。文献［9］利用Adams软件对飞机起落架的支柱进行了参数化建模，并以油孔面积为设计变量，分析了油孔面积的变化对缓冲性能的影响。文献［10］借助Optistruct结构优化平台，利用变密度法，在给出的多种工况下对飞机起落架的外筒进行了拓扑优化设计。文献［11］提出了一种施加拦阻载荷的方法，对舰载无人机拦阻着舰进行刚柔耦合多体动力学仿真分析。然而，上述学者仅仅是在给定工况条件下，对其进行了优化设计，工况条件一般是通过试验、与理论计算得来的，用这两种方法获得的工况条件不仅耗时费力，而且大都是针对特定机型而获得的工况条件，用这种工况条件进行优化设计往往不准确。

通过Adams软件建立了飞机前起落架的动力学仿真模型，然后对其进行落震仿真分析，并利用动力学分析所获得的工况条件，利用变密度法与优化准则法以扭力臂的单元相对密度作为设计变量、以柔度最小作为目标函数对飞机前起落架的扭力臂结构进行了拓扑优化设计。

2 计算方法

在以往对飞机前起落架的研究基础上，利用Adams软件建立了飞机前起落架的落震仿真模型，并对其进行了动力学分析，最后基于变密度法与优化准则法对飞机前起落架扭力臂结构进行了拓扑优化设计。

2.1 基于变密度法的拓扑优化数学模型

连续体结构拓扑优化模型的本质分为两类：一类是寻求0-1离散变量的组合优化，然后用离散变量0和1来描述优化模型，采用组合优化方法求解拓扑优化问题。对于规模较小的问题，这种组合优化方法拥有较强的全局寻优能力，但是当连续体结构拓扑优化的规模达到一定程度时，结果就截然相反了，不仅优化问题的求解效率大大降低，而且由于设计变量数量的增多，就会出现“组合爆炸”问题；

另一类方法，变密度法是将离散变量的优化问题松弛为连续变量的优化问题，从而使原本离散变量的设计模型变为连续设计变量的优化模型，这样就避免了“组合爆炸问题”，因此在拓扑优化设计中选用了变密度法。变密度法首先选定设计域与非设计域，然后将设计域进行有限元离散，以每个单元介于0、1之间的单元相对密度作为设计变量，将一个离散的优化问题转变成更易求解的连续性优化问题［12］。在运用变密度法对研究对象进行拓扑优化设计时，需要引入惩罚因子对出现的中间密度值进行惩罚，使中间密度值向0/1两端聚集，在此次优化设计中采用的材料插值模型为SIMP法，SIMP法用公式表示为：

式中：E（bi）、E1、E0、bi、p-插值之后的弹性模量、孔洞部分材料的弹性模量、实体部分材料的弹性模量、单元相对密度，取值为0的代表无材料，取值为1的代表有材料和惩罚因子。以单元相对密度作为设计变量，结构的最小柔度作为目标函数，基于变密度法SIMP法的拓扑优化数学模型可以表示为：

式中：B、C、A、S、K、f、V、V0、bmin、bmax-设计变量、结构柔度、载荷矢量、位移矢量、结构刚度矩阵、保留的体积分数、优化后的体积、结构的初始体积、设计变量的下限值和设计变量的上限值。

2.2 优化算法

在拓扑优化问题中，用来求解的常用两种算法是优化准则法与移动渐近线法。优化准则法是以目标函数和约束函数的一阶导数之商来更新单元密度变量，此法最大的特点是适用于设计变量比较多与约束条件较少的情况，且迭代次数少，计算效率高，通用性较强［13，14］。移动渐近线法是利用设计点目标函数与其一阶导数构建一个近似凸函数去逼近实际的隐函数，此法适用于求解较为平滑的非线性优化问题，缺点是寻找近似函数比较困难。此次拓扑优化问题中，选择更易求解的优化准则法。

式中：w、n、η、λ-正的移动步长、迭代次数、阻尼系数和拉格朗日乘子。

3 落震仿真模型的建立

根据某军用A型飞机前起落架的实际尺寸，建立了前起落架的三维模型，并对此模型进行了合理的简化，然后将模型保存为X_T格式并将其导入到Adams中，接着对其添加必要的运动副与力的加载。用一个小球来模拟机身的重量，并将小球与机身用固定副相连接，因为活塞杆与外筒之间既可以相互转动，又可以相互轴向移动，故在二者之间添加圆柱副；考虑到防扭臂可相对于外筒及活塞杆转动，通过旋转副分别与之连接；起落架在着陆过程中，斜撑杆处于锁死状态，故使用固定副与外筒连接；在落震仿真过程中，考虑到轮胎有一定的转速，故在轮胎与轮轴之间通过旋转副相连接；在跑道与机轮之间定义了一个接触，并用impact函数来计算接触力；在起落架缓冲系统中主要是对空气弹簧力、油液阻尼力与结构限制力的加载，由于在Adams中无法模拟油液与空气，故在外筒轴线上一点MARKER_7与活塞杆最顶端MARKER_8之间通过添加单作用力来分别模拟空气弹簧力、油液阻尼力与结构限制力。根据空气弹簧力的公式，在空气弹簧力的函数定义窗口输入其表达式：initial_air_pressure*pneumatic_area*（initial_air_volume/（initial_air_volume-pneumatic_area*（uppercylinder_length-DM（MARKER_7，MARKER_8））））**polytropic_exponent，由于油液阻尼力与活塞杆的运动方向、运动速度都有关系，因此利用VR函数与IF函数模拟了油液阻尼力，部分函数表达式为：IF（VR（MARKER_7，MARKER_8）：表达式1，0，表达式3），其中表达式1与表达式3为油液阻尼力的公式；最后对结构限制力进行模拟，从结构限制力的公式可以看出，结构限制力是一个分段函数，因此利用IF函数模拟了结构限制力，其函数表达式为：

stopper_stiffness*（uppercylinder_length-DM（MARKER_7，MARKER_8））*IF（uppercylinder_length-DM（MARKER_7，MARKER_8）：1，0，0），如图1所示飞机前起落架的落震仿真模型。

图1 前起落架的落震仿真模型Fig.1 Simulation Model of Nose Landing Gear

4 仿真分析

在飞机前起落架着陆之前需做一个静平衡，让飞机前起落架在前两秒时间里是静力学分析，两秒之后飞机前起落架才开始下落，初始仿真参数由某军用A型飞机前起落架设计参数通过计算得到［15］，部分数据如表1所示。然后设定仿真时间为4s，迭代子步为400步，最大压缩行程为600mm，如图2所示为飞机前起落架空气弹簧力随时间的变化曲线。

表1 初始仿真参数及取值Tab.1 Initial Simulation Parameters and Values

图2 空气弹簧力随时间变化规律Fig.2 Variation Law of Air Spring Force with Time

图2 中可以看出，由于前2.74s属于飞机空中降落阶段，此时空气弹簧力是一个定值约为17792 N；在2.74s时，此时飞机机轮刚好接触地面，缓冲器开始压缩，经过0.32s后空气弹簧力达到最大值约为84293 N，之后缓冲器开始回弹，空气弹簧力逐渐下降，在3.58s时达到最小值约为37600 N，完成了一次循环往返，之后空气弹簧力震荡幅度越来越小，最终趋于稳定，如图3所示为缓冲器行程随时间变化曲线。

从图3中可以看出，在前2.74s属于飞机空中降落阶段，此时缓冲器并没有开始压缩，缓冲器的行程为0，在2.74s后缓冲器开始压缩且在3.06s时缓冲器的行程达到最大值，最大值约为439 mm；随后缓冲器开始释放能量，行程开始回弹，由图可知当行程回弹到约239 mm时，完成了第一次往返行程，又开始第二次的循环往返行程时间约为0.84s。如图4所示，为油液阻尼力随时间变化曲线。

图3 缓冲器行程随时间变化规律Fig.3 Variation of Buffer Stroke with Time

图4 油液阻尼力随时间变化规律Fig.4 Oil Damping Force Changing with Time

从图4可以看出，曲线在2.74s后油液阻尼力迅速增大，与缓冲器行程、空气弹簧力增大的时刻一致，在t=2.78s时油液阻尼力达到最大值为220300 N，之后迅速减小，成正负交替震荡，最终趋于零，说明外筒与活塞杆之间相对运动速度迅速趋于减小。如图5所示，为缓冲器轴向力的变化规律。

图5 缓冲器轴向力的变化规律Fig.5 Variation Law of Buffer Axial Force

在图5中，可以看出t=2.18s时，最大轴向力为Py=123770 N，且轴向力最终趋于稳定，在拓扑优化设计中选取最大轴向力作为危险工况对扭力臂进行拓扑优化。最大起转载荷为［16］：

式中：Ay、Ax-最大起转垂直载荷和最大起转水平载荷。

由此计算出缓冲器所受扭矩为Mt=Ax*r=6280N*m，r-扭力臂与活塞杆连接处轴线到活塞杆轴线的距离。

5 拓扑优化设计

在进行拓扑优化之前，先对前起落架进行静力学分析，前起落架材料为40CrMnSiMoVA，强度极限为1760 MPa，弹性模量为200 GPa，泊松比为0.3。由于上扭力臂与外筒之间是相互铰接在一起，而外筒是固定不动的，故在上扭力臂两端施加固定约束，然后在活塞杆上施加大小为6280N*m的扭矩，通过有限元方法得到上下扭力臂的变形及应力分布情况，如图6所示。

图6 优化前扭力臂变形、应力图Fig.6 Deformation and Stress Diagram of Torsion Arm Before Optimization

从（a）图中可以看到扭力臂的最大变形量约为0.84mm，变形量很小，且从（b）图中可以看到，扭力臂的最大应力约为769MPa，最大应力值出现在扭力臂的边角部分，中部区域应力值较小，远小于材料的强度极限，这说明中间区域扭力臂的刚度冗余很大，具有很大的优化设计空间。考虑到扭力臂的优化设计空间较大，进一步对扭力臂进行拓扑优化设计，在拓扑优化设计过程中需要选择设计域与非设计域，此次优化设计过程中设计域选择为除去扭力臂吊耳部分的中间区域，体积保留分数选择为0.7，优化后的拓扑优化结果，如图7所示。

从图7中可以看到，拓扑优化结果比较清晰，扭力臂中间部分形成了明显的开槽区域，按照拓扑优化生成的结果，在三维软件中，重新对起落架进行精确建模。

图7 拓扑优化结果云图Fig.7 Cloud Chart of Topology Optimization Results

接下来对重新建立后的模型，进行强度设计校核，利用有限元方法得到了扭力臂优化后的变形与应力分布结果，如图8所示。

图8 优化后扭力臂变形、应力图Fig.8 Deformation and Stress Diagram of Torsion Arm after Optimization

从图8中可以看到优化后的扭力臂最大变形量约为0.88mm，比优化前的变形量仅增加了0.04mm；优化后的最大应力值约为842MPa，比优化前的应力值增加了73MPa，增加的百分比约为9.4%，且优化后的应力值仍在材料的强度极限范围内，故优化后的结构满足强度要求；优化后扭力臂的质量为7.96kg，优化前扭力臂的质量为10.68kg，优化后扭力臂的质量比优化前的质量减少了2.72kg，减少的百分比为25.4%，轻量化效果明显。

6 结论

利用Adams软件建立了飞机前起落架的落震仿真模型，并进行仿真分析，得到了飞机前起落架缓冲系统轴向载荷及行程的变化规律，最后在动力学分析的基础上基于变密度法与优化准则法对飞机起落架进行了拓扑优化设计，得到如下结果：

（1）在降落阶段，此时空气弹簧力是一个定值，当飞机着陆后缓冲器开始压缩，空气弹簧力迅速增大，增大到峰值后，缓冲器开始回弹，空气弹簧力迅速下降，之后缓冲器重复此过程，空气弹簧力来回震荡，且震荡幅度越来越小，最终趋于稳定；缓冲器行程在空中降落阶段为0，此后变化规律与空气弹簧力变化规律一致；在空中降落阶段，油液阻尼力也为0，在之后的0.3s内，油液阻尼力震荡幅度很大，此后震荡幅度很小，直至趋于稳定；

（2）在动力学分析的基础上，基于变密度法与优化准则法以单元相对密度作为设计变量，以结构柔度最小作为目标函数，在满足强度与刚度的条件下，对扭力臂结构进行了拓扑优化设计，得到的拓扑优化结果清晰明了且易于加工，优化后扭力臂质量减小了2.72kg，减少的百分比为25.4%，轻量化设计效果明显。