垂直着陆中直升机旋翼动力学行为研究

胡国才，吴靖，刘湘一，刘书岩

海军航空大学 航空基础学院，烟台 264001

目前大、中型铰接式或弹性轴承旋翼直升机广泛采用刚度低、阻尼大的液压阻尼器，提高旋翼摆振面阻尼以避免直升机发生地面共振[1]。液压阻尼器的阻尼力与活塞速度或桨叶摆振速度的平方近似呈比例关系，为了限制使用中阻尼器的动态载荷，以满足阻尼器和旋翼结构的强度及寿命要求，阻尼器内部一般设置定压安全活门，当阻尼力超过活门开启压力时安全活门打开，从而限制载荷的进一步增加[2-3]。设置安全活门带来的主要问题是直升机起降及地面运转时，若旋翼发生大的扰动安全活门开启后，阻尼器的等效阻尼会迅速下降，对地面共振产生不利影响[4]。因此，在正常使用中不会发生地面共振的直升机可能在粗暴着陆时因过大的初始扰动而发生地面共振，这在实际使用中仍时有出现。

直升机地面共振详细设计要求GJB720.5—89规定，如果直升机采用显著影响地面共振的非线性元件，则应给出不发生地面共振的初始扰动幅值，并在设计使用包线内可能遇到的扰动下应不发生地面共振。带有安全活门的液压阻尼器具有严重的非线性，采用等效线性系统，可方便地确定满足地面共振稳定性的初始扰动幅值。但在设计使用包线内可能遇到的扰动幅值，在直升机设计时难以预估，这为直升机地面共振设计带来了不确定性。虽然也可以通过半主动或主动控制的方法使直升机工作在全稳定状态，如直接使用半主动式阻尼器磁流变阻尼器[5]，或者增加机体状态反馈[6]等来达到对直升机地面共振的自适应控制[7]。但是，目前使用最多的还是结构相对简单的被动式阻尼器，这就要求尽量提高阻尼器的等效阻尼避免直升机产生地面共振，但要满足结构强度与寿命的要求，必然付出结构重量代价。要合理解决以上两者之间的矛盾，前提是需要准确预估直升机使用包线内可能遇到的扰动大小。

经验表明，直升机单轮粗暴着陆或着舰产生的扰动可能最为严重。现有关于直升机着陆和着舰方面的研究，重点关注起落架及机身结构的载荷、强度及疲劳寿命。许多学者采用多体动力学模型，用经典方法预估直升机着陆载荷[8]；用刚-柔混合有限元模型计算着陆载荷及机体和起落架组件的应力[9-10]。文献[11]在全局有限元模型中改进了关键部位模型，研究粗暴着陆时关键部位的疲劳应力。但关于直升机着陆或着舰时激起的旋翼扰动及动力稳定性问题的研究报道很少。

本文通过直升机着陆动力学建模与仿真，从理论上分析垂直着陆时起落架撞击地面对旋翼动力学行为的影响，希望为预估粗暴着陆时阻尼器的扰动幅值提供一种有效的方法。

1 分析模型

假设直升机限定在垂直着陆状态，考虑质心垂向运动、机体绕质心的滚转和俯仰等3个自由度，不考虑质心在水平面内的移动和机体的偏航自由度。因机体刚度比起落架刚度要大得多，可将机体看做刚体。对于铰接式或弹性轴承旋翼来说，对阻尼器轴向运动及地面共振起主要作用的是桨叶绕挥舞铰和摆振铰的转动，故本文只考虑刚硬桨叶的挥舞、摆振和变距自由度。起落架轮胎底面接触地面后，轮胎和缓冲支柱压缩产生弹性力和阻尼力并作用到机体上，对机体产生力和力矩，引起机体和旋翼的扰动运动。直升机垂直着陆动力学分析模型如图1所示。

图1 直升机垂直着陆模型Fig.1 Vertical landing model of helicopter

图1中，Ogxgygzg为地面坐标系，也是惯性坐标系；Ofxfyfzf为机体坐标系，跟随机体运动，原点Of位于直升机质心处，Ofxf指向机身尾部；Ohxhyhzh为桨毂不旋转坐标系，跟随机体运动，原点Oh位于桨毂中心，Ohxh指向机身尾部；Obxbybzb为桨叶坐标系，跟随桨叶运动，原点Ob位于弹性轴承处，Obxb固定在桨叶上。设旋翼桨叶片数为Nb，旋翼转速为Ω(俯视逆时针旋转)，第k片桨叶的方位角为ψk、挥舞角为βk(向上为正)、摆振角为ζk(逆旋转方向为正)、变距角为θk，弹性轴承的外伸量为e，机体滚转和俯仰欧拉角为γ和ϑ，机体滚转和俯仰角速度为ωx和ωy，旋翼轴到质心距离为xc，桨毂中心到质心高度为h。

任一片桨叶上距离弹性轴承r处的叶素在机体坐标系中的矢径为

ρ=[xf,yf,zf]if

(1)

式中：if为机体坐标系单位向量，且

(2)

叶素在地面坐标系中的矢径为

ρ=[xg,yg,zg]ig

(3)

式中：ig为地面坐标系单位向量，且

(4)

叶素在地面坐标系中的速度和加速度可由式(3)分别对时间求一次和二次导数得到。把叶素在惯性坐标系中的速度及加速度转换到桨叶坐标系中，可以计算叶素的气动力及惯性力。

着陆时起落架撞击地面会引起机体和旋翼的扰动运动，为了提高分析精度，应计入非定常气动力影响。对于旋翼低频扰动运动，Pitt/Peters动力入流模型可反映非定常气动力的影响[12]，这在旋翼/机体耦合动稳定性研究中得到了验证[13]。考虑到下降速度对旋翼诱导速度的影响，采用扩展的动力入流模型[14]。设桨盘平面上任一点的无因次诱导速度为

v=v0+vs(r+e)sinψ+vc(r+e)cosψ

(5)

式中：v0为平均诱导速度；vs和vc分别为横向和纵向入流系数。Pitt/Peters动力入流方程为

(6)

式中：M为显在空气质量矩阵；L为入流增益矩阵；CT、CL及CM分别为旋翼总的气动升力系数、对桨毂中心的气动滚转及俯仰力矩系数。

前飞状态时矩阵L为

(7)

式中：

直升机以速度vz垂直下降时，前进比μ=0，入流角as=-90°，上述系数改写为：vT=v0-vz，vm=2v0-vz，α=90°。那么垂直下降时的矩阵L为

(8)

直升机接近地面时产生的地面效应，会降低旋翼的诱导速度，增加旋翼升力[15]。旋翼总距一定的情况下，旋翼升力增加会增大旋翼平均挥舞角和摆振角，但对旋翼周期挥舞和周期摆振的影响相对较小，因此可不考虑地面效应对旋翼扰动运动的影响。

采用弹性轴承旋翼系统，液压阻尼器安装时一般将一端用球铰与桨毂连接，另一端用球铰与桨叶或过渡件连接，这样桨叶绕弹性轴承的挥舞、摆振和变距运动都会影响到阻尼器速度，故计算阻尼器速度时需计入几何耦合的影响[16]。液压阻尼器的力-速度曲线由试验数据得到，但为了仿真计算的方便，一般采用双折线来拟合试验数据获得理论计算公式。

起落架触地时机体姿态角不大的情况下，沿缓冲支柱的轴向载荷对机体起主要作用，本文暂不考虑轮胎触地后产生的侧向力和纵向力。

根据达朗贝尔原理导出直升机垂直着陆的动力学方程，以一阶向量表示为

(9)

机体的运动学方程为

(10)

图2 直升机起落架模型Fig.2 Model of helicopter landing gear

为便于计算起落架着陆载荷，建立起落架参考坐标系。起落架触地后压缩或伸展时，缓冲支柱和轮胎均会产生非线性弹性力和阻尼力，分别以非线性弹簧S1、S2及非线性黏壶D1、D2表示，起落架模型如图2所示。前起落架用下标“N”表示，左侧主起落架用下标“ML”表示，右侧主起落架用下标“MR”表示。

起落架参考点Or与直升机质心处于同一铅垂线上，参考点垂直位移与直升机质心垂直位移之间的关系为zr=-z。起落架的垂直位移取决于参考点的垂直位移zr和机体姿态角。假定参考点到前起落架的距离为xN，到主起落架连线的距离为xM，主起落架轮距为2yM。

则起落架的垂向位移为

(11)

起落架的垂向位移对时间求导，可得起落架的垂向速度。

(12)

式中：Nw为单个起落架的轮胎数量。

缓冲支柱和轮胎的弹性力取决于压缩量，为提高预估精度，采用静压缩试验数据。缓冲支柱的阻力包含油液流经阻尼孔产生的阻力fd1h和活塞摩擦力fd1c。设ρ为油液密度、Cd为流量系数、A0为节流阀座的承压面积、As为压缩行程的阻尼孔面积、Ar为伸展行程的阻尼孔面积，则流体阻力为

(13)

活塞摩擦力fd1c大小与活塞的摩擦系数及缓冲支柱内部空气压力有关，其方向与活塞运动方向相反。影响摩擦力的因素较多，需要由试验确定。由于缺少摩擦力试验数据，本文暂不考虑摩擦力的作用。轮胎的阻尼力采用经验公式[17]：

(14)

式中：Cw为轮胎当量阻尼系数。

起落架作用在机体上的垂向力Fz、滚转力矩Mx和俯仰力矩My，即各起落架缓冲支柱的弹性力和阻尼力之和。仅考虑起落架的垂向载荷时，作用在机体上的力和力矩分别为

(15)

起落架过载系数定义为

ngi=(fs1i+fd1i)/p0ii=N,ML,MR

(16)

式中：p0i为直升机在零升力停机状态时各个起落架的地面支反力。

2 着陆仿真与结果分析

2.1 计算方法

采用4阶Runge-Kutta法对直升机着陆运动方程进行数值积分，模拟直升机垂直着陆过程中旋翼桨叶、机体及起落架的动态响应。给定不同的高度和旋翼总距，可获得不同的着陆速度；设置不同的周期变距可使直升机产生不同的姿态角。首先将直升机高度固定，设置旋翼总距和周期变距，对运动方程进行初值积分，经过一段时间后旋翼运动和机体姿态趋于平衡状态。直升机趋于平衡后，放松质心垂向自由度，使直升机垂直下降，当起落架触地时，机体在起落架支反力作用下产生扰动运动，同时激起旋翼桨叶的挥舞和摆振运动。对全部桨叶挥舞和摆振响应进行多桨叶坐标转换，得到旋翼整体模态响应，对其进行快速傅里叶变换(FFT)及滤波运算可以辨识各个模态。

旋翼和机体运动方程的积分变量为方位角ψ，起落架动力学方程的积分变量为时间t，两者之间的关系为ψ=Ωt。最终的载荷及动态响应统一以时间t或方位角ψ为自变量的形式输出。

2.2 结果分析

以某弹性轴承旋翼直升机为例，旋翼和机体的主要参数如表1所示。

不稳定流场、或驾驶员操纵不当等因素都会导致直升机单轮粗暴着陆，表现为直升机以一定倾斜角和较大速度接地。作为算例，假定让右起落架最先接地，接地速度大约为2 m/s。为模拟这种着陆状态，直升机需从较高的高度开始下降，驾驶杆需向右、向后拉杆。设直升机从高度为1 m(轮胎底面距离地面的平均高度)开始垂直下降着陆，旋翼总距θ0=6°、周期变距θc=-2°及θs=0.5°，重心位于中位xc=0。图3为机体滚转和俯仰角响应，图4为起落架过载系数随时间的变化曲线。

图3显示，开始阶段机体呈现向右倾斜约3.5°、抬头1°的姿态，该姿态角与右压杆(θc=-2°)及后拉杆(θs=0.5°)的配平特性相吻合。运动方程积分时，旋翼周期挥舞和周期摆振的初值设置为零，故要经过一定时间才能达到配平状态，在此过程中机体姿态角随时间会存在一定波动，机体滚转角的波动更为明显。

表1 旋翼及机体主要参数Table 1 Main parameters of rotor and body

图3 机体滚转和俯仰响应Fig.3 Body roll and pitch responses

从t=1.57 s开始下降，大约在t=2.34 s时右起落架轮胎触地，随后起落架载荷迅速增加(见图4)，机体在该载荷作用下产生左滚和低头力矩，使机体迅速向左滚转并低头，如图3所示。大约经过0.11 s，左起落架和前起落架先后(几乎同时)触地。右起落架最大过载系数为1.44，左、前起落架最大过载系数分别为1.58和1.53，比右起落架分别高9.7%和6.2%。着陆后期，机体略有向右倾斜和抬头的姿态，左、右起落架载荷也不相同，主要是因为直升机着陆后旋翼桨距没有回零的缘故。

图4 起落架过载系数随时间变化曲线Fig.4 Curves of overload factor of landing gears vs time

图5 桨叶摆振角随方位角的变化曲线Fig.5 Curves of blade lag angle vs azimuth

图5显示了各片桨叶的摆振响应随第1片桨叶方位角的变化曲线。可以看到桨叶的平均摆振角为6°左右，定常周期摆振的幅度约为0.18°。着陆过程中各片桨叶产生了不同程度的摆振扰动，扰动幅度大约在0.5°～1.2°之间。

从仿真结果看到，右起落架触地前，各片桨叶按照固定的相位差(60°)以转速Ω作定常摆振运动。一个有意思的现象是，全部起落架触地引起的扰动摆振响应与定常摆振响应叠加后，至第1片桨叶处于330°方位角时，各片桨叶均近似处于平均摆振角位置。当然，这个现象与各起落架触地时机、机体姿态等诸多因素有关，不具有普遍意义，但为下面的力学分析提供了便利。

图6表示第1片桨叶处于330°方位(坐标轴xf为0°、方位角按顺时针确定)，前、左起落架接地后对机体产生向后及向右的力矩，机体扰动使桨毂中心产生向后及向右的加速度，对桨叶施加向前和向左的惯性力作用。

从图6可以看出，此时1号和2号桨叶在惯性力矩作用下将逆旋转方向摆振，摆振角将增加(因规定了逆旋转摆振为正)，综合两个惯性力的作用，1号桨叶应比2号桨叶向后摆振更剧烈；4号和5号桨叶将顺旋转方向摆振，摆振角将减小，并且4号桨叶应比5号桨叶向前摆振更剧烈。对于3号和6号桨叶来说，由滚转和俯仰运动产生的惯性力矩起相互抵消的作用，摆振幅值小，而摆振方向与两者的大小有关，图5显示的响应曲线与上述分析结果吻合。随后的桨叶摆振运动因与机体运动、桨叶旋转及定常摆振运动等因素存在复杂关系，难以一一分析清楚。

图6 机体扰动引起的桨叶惯性力Fig.6 Blade inertial force caused by fuselage perturbation

全部起落架从2.34 s至3.0 s完成第1次压缩-回弹运动，持续时间大约为0.66 s (见图4)，之后起落架和机体处于小幅震荡。因此，在着陆扰动期间难以将摆振后退型的瞬态响应分离出来，因而也就不能根据包络线来识别其模态阻尼[18]，但可知着陆撞击激起的摆振扰动大小。

对旋翼摆振余弦分量进行2 Hz低通、4 Hz高通和3～4 Hz带通滤波可知，除了摆振后退型模态响应外，还存在其他频率的响应成分，具体如图8所示。

图7 旋翼摆振的余弦分量及2～3 Hz带通滤波后的响应Fig.7 Cosine component and 2-3 Hz band pass filtered rotor lead-lag response

很明显，幅度最大、频率最低的响应成分是由起落架触地后压缩-回弹引起的(实线)；高通滤波后仅存在频率大约为4.5 Hz的响应成分(虚线)，因摆振前进型模态频率ωζ P=4.46 Hz，可以确定这是摆振前进型响应，其幅值约为0.17°；频率为3.5 Hz持续振动的响应成分是定常周期摆振，着陆时其幅值也稍有增加(点线)。

从以上仿真结果及分析可知，单轮着陆会激起旋翼摆振后退型响应，虽然其幅值并不大，但由于存在各种频率成分的背景振动，尤其是低频、大幅值的冲击响应，液压阻尼器将处于大速度状态而打开定压活门，进而使其等效阻尼下降。图9为各片桨叶的液压阻尼器轴向速度响应，图10为该阻尼器的等效阻尼随速度幅值的变化曲线[19]，定压活门的开启速度约为0.01 m/s。

由图9看到，定常状态下阻尼器的速度大约为0.012 m/s，着陆时各阻尼器速度峰值高达0.035～0.05 m/s，已远超定压活门的开启速度，由图10看到，阻尼器的等效阻尼将随着速度的增加而急剧下降。

图8 旋翼摆振响应中的各种成分Fig.8 Components in rotor lead-lag response

图9 各阻尼器的轴向速度响应Fig.9 Axial velocity response of each damper

图10 等效阻尼随轴向速度的变化曲线Fig.10 Curves of equivalent damping vs axial velocity

3 垂直着陆对载荷及旋翼动力学影响

直升机着陆具有一定随机性，其构型、着陆场地、大气环境以及驾驶员操纵等因素都会影响着陆状态[20-21]。为模拟着陆随机性，针对本文的直升机构型，给定高度范围为0.3～1.5 m、间隔为0.3 m，在一定范围内随机给出开始下降时间、旋翼总距和周期变距，开始下降时间范围为1～2 s、总距范围为6°～7°、周期变距范围为-2°～2°。在每个高度，随机着陆5次，输出模拟参数为各起落架着陆速度和最大过载系数、桨叶挥舞和摆振动幅值以及阻尼器速度峰值，以散点图的形式显示。起落架触地瞬间的机体姿态角及角速度等动态参数，取决于着陆高度、开始下降时间及旋翼变距等因素，仿真结果如图11所示。

图11(a)和图11(b)显示了直升机不同高度随机着陆时起落架的着陆速度及过载系数，反映了着陆的严重程度。

从图11(a)和图11(b)可以看到，着陆速度及过载与着陆高度具有一定的相关性。说明着陆高度对起落架着陆速度和过载有重要影响。同时也看到，在同一高度随机着陆时，着陆速度和过载系数具有相当大的分散性，表明除了着陆高度外，触地时机体的姿态角及角速度也是影响着陆速度和过载的重要因素。主起落架着陆速度的分散性大约为0.8 m/s，基本上不随着陆高度而变。而过载系数的分散性却随着陆高度的增加而有所扩大。

图11(c)和图11(d)分别表示着陆时桨叶挥舞角和摆振角的最大动幅值。

图11 仿真结果图Fig.11 Results of simulation

从图11(c)和图11(d)可见，随着陆高度的增加，挥舞角和摆振角的最大动幅值的分散性明显增大。说明着陆高度越大，随机着陆时桨叶可能会遇到更大的扰动。与挥舞相比，桨叶摆振扰动的分散性更为显著，原因是摆振面气动力远小于挥舞面气动力，因此着陆时机体扰动引起的摆振气动阻尼也相应地小得多，而扰动产生的惯性力对桨叶摆振运动的作用要大于对挥舞运动的作用。从图11(d)看到，高度为0.3 m随机着陆时，摆振角幅度在0.41°～0.61°之间，而在1.5 m高度随机着陆时，摆振幅度范围扩大至0.67°～1.7°。

图12表示随机着陆时液压阻尼器的轴向速度峰值随着陆高度的变化。

图12 液压阻尼器轴向速度峰值随着陆高度的变化Fig.12 Maximum of hydraulic damper axial velocity vs landing altitude

由图12看到，在不同高度随机着陆时，阻尼器轴向速度峰值也存在很大的分散性。由于阻尼器与桨叶挥舞、摆振及变距自由度存在几何耦合，其轴向速度主要取决于桨叶的摆振速度外，还与桨叶挥舞和变距运动有关。因此阻尼器的轴向速度峰值与图14所示的桨叶最大摆振动幅值之间并不是严格的对应关系。不同高度随机着陆时，阻尼器轴向速度峰值的最小值处于24～30 mm/s范围内，最大值处于45～74 mm/s范围内，说明粗暴着陆时，阻尼器的轴向速度峰值会显著增加。

4 结 论

通过建立直升机垂直着陆动力学模型并进行数值仿真，研究了直升机垂直着陆中起落架、机体及旋翼的动力学行为，分析了着陆状态对动态响应参数的影响规律，得到以下主要结论：

1) 直升机着陆时，起落架撞击地面引发的机体滚转和俯仰扰动是激起旋翼动态响应的主要原因。

2) 起落架触地后的压缩-回弹周期中，将激起各片桨叶不同的摆振幅值，并形成旋翼摆振后退型响应。

3) 随机着陆时，起落架着陆速度及过载系数与着陆高度具有明显的相关性，而机体姿态角及角速度的不同使得同一高度着陆时起落架的着陆速度及过载具有一定分散性。

4) 随着陆高度增加，着陆姿态对桨叶挥舞和摆振动幅值的影响增加，尤其对桨叶摆振动幅值的影响更为显著，随着着陆高度增加，摆振动幅值的分散性也增加。

5) 不同高度随机着陆时，阻尼器轴向速度峰值呈现很大的分散性，且由于阻尼器与桨叶挥舞、摆振及变距自由度存在几何耦合，使得其轴向速度峰值的分散性随着着陆高度的变化规律不明显。