水平着陆飞行器着陆架系统耦合压缩触地过程仿真分析

刘 刚，郑宏涛，李 洋

(中国运载火箭技术研究院研究发展中心， 北京 100076)

0 引言

水平着陆飞行器着陆架触地是一个短时间、高动态的复杂运动过程。为了实现飞行器安全着陆和触地，保证飞行器在此过程中不后翻、尾部不擦地、轮胎和减震器压缩量不超限，需要对触地过程进行精确的仿真分析，才能确定或优化相关的设计参数。着陆架触地过程中轮胎和减震器的压缩运动存在耦合因素，因此轮胎压缩量及其受力的求解成为问题的关键环节。高泽迥[1]给出了着陆架减震器和轮胎的力学模型，但是未给出压缩量的求解方法。罗琳胤等[2]采用商业软件对起落架多体动力学进行建模与分析，但是软件中的压缩量求解方法未知。陈丽城等[3]采用Matlab/Simulink实现了无人机地面动力学建模及分析，但是其模型中不包含减震器和轮胎压缩量的耦合因素。针对水平着陆飞行器着陆架系统耦合压缩触地过程仿真分析问题，本文采用一种基于瞬时平衡条件的求解方法，先求解减震器压缩量，再求解轮胎压缩量，进而求得轮胎法向力，实现水平着陆飞行器触地过程的仿真分析。

1 着陆架系统耦合压缩模型

飞行器摇臂式着陆架系统示意图如图1所示。其中AB为悬臂梁，与飞行器机体固连。ACD为摇臂，在A与悬臂梁铰接，D为轮轴。BC为减震器，在B与悬臂梁铰接，在C与摇臂铰接。

图1 摇臂式着陆架系统示意图Fig.1 Schematic figure of a rock arm undercarriage system

fa=fa(S)

(1)

(2)

(3)

轮胎法向力N是轮胎压缩量δ的函数,滚动摩擦力f与法向力成线性关系，如式(4)、式(5)所示，式中cf为轮胎的滚动摩擦系数。

N=N(δ)

(4)

f=cfN

(5)

在轮胎下边沿与着陆场地面接触后，飞行器开始触地过程。由于飞行器在触地时具有一定的垂直方向的下沉速度，会造成轮胎迅速压缩，轮胎压缩后轮轴受到支持力和滚动摩擦力。这两个力对摇臂产生一个对A点的逆时针力矩，当这个力矩增大到一定程度时，摇臂会逆时针转动。摇臂的逆时针转动一方面会造成减震器压缩；另一方面会改变着陆架系统的几何构型，使轮轴D的高度增大。轮轴D的高度增大会造成轮胎压缩量减小。因此着陆架的减震器和轮胎的压缩量是耦合的。

2 基于瞬时平衡条件考虑着陆架系统耦合压缩的触地过程仿真分析方法

2.1 瞬时平衡条件推导

摇臂的运动过程满足绕A点转动的微分方程：

(6)

其中，JACD为摇臂的转动惯量，ωACD为其转动角速度，∑MA为摇臂受到的绕A点的力矩之和。

通过直接求解这个微分方程对着陆架系统进行仿真具有一定的困难。一方面，压缩量的耦合因素使得每一步仿真中的受力和力矩难以确定；另一方面，如果直接求解微分方程，由第k-1步递推到第k步时，可能出现减震器受摩擦力和阻尼力与其运动趋势不符合的情况，即使仿真步长比飞行器六自由度仿真步长低一个数量级，求解结果也会发生振荡或跳变。为解决此问题，本文采用一种基于瞬时平衡条件的着陆架系统耦合压缩量和受力的仿真分析方法。

对上述摇臂的运动方程进行分析。一般情况下，摇臂的转动惯量与角加速度的乘积相对于来自轮胎和减震器的力矩是小量，其数值一般低3个以上数量级。因此可忽略此小量，即得到摇臂的瞬时平衡公式：

∑MA=0

(7)

2.2 着陆架减震器运动趋势判断方法

摇臂所受的对A点的力矩包括5部分：轮胎法向力N对应的力矩MN，轮胎滚动摩擦力f对应的力矩Mf，减震器弹簧力fa对应的力矩Mfa，减震器摩擦力ff对应的力矩Mff，减震器阻尼力fd对应的力矩Mfd。由于减震器摩擦力和阻尼力的方向与减震器运动方向相反，因此需要先判断减震器的运动方向，再分不同的情况以不同的方法计算合力矩。

在第k步仿真中，来自飞行器六自由度运动方程的相关参数(包括飞行器质心高度hk，俯仰角φk，滚转角φk，质心在飞行器的坐标Xgk、Ygk、Zgk)被更新。此时，假设减震器压缩量暂时维持与第k-1步的值不变，则根据第k-1步仿真中轮轴D的坐标Xdk-1、Ydk-1、Zdk-1和减震器压缩量Sk-1以及相关几何关系可计算出一个轮胎压缩量δtemp：

δtemp=δtemp(hk,φk,φk,Xgk,Ygk,Zgk,

Xdk-1,Ydk-1,Zdk-1,Sk-1)

(8)

3)减震器压缩量不变情况：若不满足上述两种情况，则减震器压缩量维持第k-1步的值不变。

2.3 减震器和轮胎的压缩量耦合求解方法

分别针对上述3种情况，给出第k步仿真中减震器和轮胎的压缩量的耦合求解方法。

(1)减震器压缩量不变情况

此种情况最简单，压缩量按式(9)、式(10)计算：

Sk=Sk-1

(9)

δk=δtemp

(10)

(2)减震器压缩量增大情况

(11)

其中,step为仿真步长，与飞行器六自由度仿真步长相同。

如图1所示，减震器压缩量改变后，着陆架系统的几何关系将发生改变，C点变为C1，D点变为D1。因此需要根据S*k的值重新计算几何关系，计算出轮轴D1的坐标Xd*k、Yd*k、Zd*k，重新计算轮胎压缩量δ*k：

δ*k=δ*k(hk,φk,φk,Xgk,Ygk,Zgk,

Xd*k,Yd*k,Zd*k,S*k)

(12)

对A点的和力矩为：

∑M*A=M*N+M*f-(M*fa+M*ff+M*fd)

(13)

注意到，每对应一个假定的减震器压缩量ΔS*,可得到一个∑M*A=∑M*A(ΔS*)，这是一个以ΔS*为自变量的一元函数。按前文推导的瞬时平衡条件可得：

∑M*A(ΔS*)=0

(14)

式(14)是一个以ΔS*为自变量的一元非线性方程，可用割线法[4]等数值法求解，解方程时需要注意限制自变量的取值在合理范围。

其解即为第k步的减震器压缩量增量ΔSk，则第k步减震器压缩量为：

Sk=Sk-1+ΔSk

(15)

根据Sk重新更新几何关系，包括轮轴的坐标Xdk、Ydk、Zdk，各杆的长度和各角的角度等，进而可计算出第k步的轮胎压缩量δk。

δk=δk(hk,φk,φk,Xgk,Ygk,Zgk,

Xdk,Ydk,Zdk,Sk)

(16)

(3)减震器压缩量减小情况

此种情况求解方法与减震器压缩量增大情况类似，唯一不同的是求解的瞬时平衡方程如式(17)：

∑MA(ΔS) =Mfa-(Mff+Mfd+MN+Mf)=0

(17)

2.4 飞行器触地过程仿真分析方法

根据2.3节的方法分3种情况可求得第k步的减震器压缩量Sk和轮胎压缩量δk。根据轮胎法向力与轮胎压缩量的函数关系或插值表可求得轮胎法向力Nk。由滚动摩擦系数可求解滚动摩擦力fk。在仿真第k步，把飞行器视为一个整体，将Nk、fk及其对飞行器质心的力矩加入飞行器六自由度运动方程并进行求解，即可实现考虑着陆架系统耦合压缩的触地过程仿真分析。

3 考虑着陆架系统耦合压缩的触地过程仿真分析算例

采用本文给出的方法，选取触地下沉速度-1.51m/s、触地俯仰角11.2°、触地滚转角-0.02°的典型工况，对水平着陆飞行器的触地过程进行仿真分析，仿真步长取1ms。仿真结果如图2～图9。

图2 左减震器压缩量Fig.2 Compressing amount of the left shock absorber

图3 右减震器压缩量Fig.3 Compressing amount of the right shock absorber

图4 左主轮压缩量Fig.4 Compressing amount of the left main tire

图5 右主轮压缩量Fig.5 Compressing amount of the right main tire

图6 飞行器俯仰角Fig.6 Pitch angle of the vehicle

图9 飞行器尾尖与地面距离Fig.9 Distance between the vehicle tail and ground

由如图2～图9可见，左、右主轮都被弹起一次，并进行了二次触地。减震器压缩量和轮胎压缩量都经历了从零到最大再回落并稳定的过程。飞行器的俯仰角、滚转角和下沉率都逐渐收敛于零附近，并保持稳定。在触地过程中飞行器未发生后翻、尾部擦地、轮胎或减震器压缩量超限现象。

4 结论

针对水平着陆飞行器触地过程六自由度运动与着陆架减震器和轮胎的耦合压缩运动联合仿真问题，将复杂的着陆架的机构多体动力学进行简化，推导了摇臂式着陆架的力矩瞬时平衡条件，给出了着陆架减震器运动趋势判断方法，通过求解以减震器压缩量增量为自变量的一元非线性方程求得减震器压缩量，进而求得轮胎压缩量和轮胎受力，实现飞行器触地过程仿真分析。该方法可在飞行器六自由度仿真模型基础上直接扩展，无需改变其仿真步长，着陆架模型求解无需依赖商业软件，可对水平着陆飞行器触地过程相关参数和性能进行有效的仿真评估和验证。

参考文献

[1] 高泽迥.飞机设计手册第14册：起飞着陆系统设计[M].北京：航空工业出版社，2002.

[2] 罗琳胤，边宝龙.飞机起落架缓冲性能仿真分析[J].机械设计，2012,29(4)：56-58+62.

[3] 陈丽城, 李春涛, 张孝伟, 等. 无人机地面动力学建模及分析[J]. 计算机仿真, 2016 (6): 13-18.

[4] 颜庆津.数值分析[M].北京：北京航空航天大学出版社，2000.