董 倩， 张献民， 包伊婷， 程少锋， 张宇辉
(1. 中国民航大学 交通科学与工程学院， 天津 300300；
2. 天津大学 建筑工程学院， 天津 300057；
3. 南京航空航天大学 民航学院， 南京 210016)
准确确定飞机滑跑过程中对跑道的动荷载是研究道面结构动力学特性的前提。根据研究目的不同，国内外将飞机荷载简化为静载[1-3]、移动恒载[4-8]、移动简谐荷载[9]及冲击荷载[10-11]，并在此基础上开展了有益的研究。然而，飞机滑行过程中的动荷载除受自身载重外，主要受到道面平整度及飞机升力的影响。因道面平整度是具有零均值、各态历经平稳的高斯随机过程[12]，飞机升力亦随滑行速度变化明显，故飞机在起降滑行过程中产生的动荷载随时间与空间位置不断变化，亦具备随机性。因此，有必要针对飞机随机动荷载的时域及频域特性展开研究。
国内外围绕飞机起降滑行过程中的随机动荷载开展了一些研究。魏保立等[13]建立二自由度飞机-道面结构的耦合振动分析模型，计算了飞机随机动荷载，分析了荷载作用下道面的随机振动响应，结果表明，飞机随机荷载中的高频作用力更容易引起道面结构的振动；
孙璐等[14]在描述机场道面平整度随机特性的基础上，将飞机简化为二自由度振动系统，利用频响函数获得了飞机对道面的随机动压力谱，并证明了随机动压力同样具备正态、零均值分布特征；
Li等[15]为了克服离散傅里叶变换在分析道面平整度过程中无时域功能的局限性，提出了利用小波变换分析道面平整度与飞机动荷载的方法。道路工程中，孙吉书等[16]建立了四轴重载车辆七自由度模型，分析了车-路的耦合作用，编写Matlab程序计算了车辆动荷载；
李倩[17]基于七自由度车辆模型研究了车辆动荷载的特性，获取了不同路面平整度等级、不同速度、不同载重下的车辆随机动荷载；
Cole等[18]对英格兰的各级公路进行实地测量，结合对常见运输车辆的振动情况统计，发现车辆动载频率在1～16 Hz。
上述研究为本文的研究提供了借鉴，但围绕飞机四自由度模型、飞机随机动荷载的频域特征及影响因素等方面的研究有待进一步开展。基于此，本文通过建立飞机主起落架四自由度模型、分析飞机-道面的耦合作用，求解了升力影响下飞机的随机动荷载，通过求解飞机动荷载的功率谱密度研究了其频域特性，并探讨了影响飞机动荷载功率谱密度的主要因素。研究结果旨在为飞机设计参数优化、机场道面设计及运维奠定基础。
常用民航客机前起落架承担的荷载通常不超过飞机整机荷载的10%，仅建立飞机主起落架模型，如图1所示。对图1所示的飞机模型做如下假设：① 飞机主起落架模型的重心在主起落架中心位置，飞机被简化为由悬挂质量、非悬挂质量、弹簧、阻尼组成的质阻弹振动系统，悬挂质量与非悬挂质量之间由弹簧连接；
② 飞机具有四个自由度：机身质心处垂直位移、机身侧倾转角以及左右主起落架位置处的垂向运动；
③ 悬挂部分的刚度、阻尼和非悬挂部分的刚度、阻尼均为常数；
④ 轮胎与道面始终处于接触状态，道面平整度作为飞机系统输入响应。
图1 飞机主起落架四自由度模型Fig.1 Model of aircraft main landing gear with four freedom degrees
根据达朗贝尔原理，机身质量M0垂向运动的振动微分方程见式(1)。其他三个运动的振动微分方程见式(2)～式(4)。
(1)
机身质量M0侧倾运动
(2)
左主起落架非悬挂部分垂向运动
(3)
右主起落架非悬挂部分垂向运动
(4)
式中：M0为悬挂质量，主要包括机身质量、机翼与缓冲器外筒质量(kg)；
Z0为飞机质心处的垂向位移(m)；
KXHZ为悬挂部分刚度(N·m-1)；
CXHZ为悬挂部分阻尼(N·s·m-1)；
KLHZ为非悬挂部分刚度(N·m-1)；
CLHZ为非悬挂部分阻尼(N·s·m-1)；
Z3为非悬挂部分垂向位移(左) (m)；
Z4为非悬挂部分垂向位移(右) (m)；
Z5为轮胎绝对垂直位移 (左) (m)；
Z6为轮胎绝对垂直位移 (右) (m)；
ly为后轴左右轮距的一半(m)；
Jy为飞机侧倾转动惯量(kg·m2)；
θy为飞机侧倾角位移；
Y为飞机升力(N)；
M1为非悬挂部分质量(左) (kg)；
M2为非悬挂部分质量(右) (kg)，M1=M2，主要包含刹车装置、轮胎等质量。
飞机左右两个主起落架位置处轮胎的受力见式(5)与式(6)
(5)
(6)
式中：Fd1为左飞机轮胎处的动荷载(N)；
Fd2为右飞机轮胎处的动荷载(N)；
M1、M2含义与式(3)、(4)中相同。
飞机在滑跑过程中机翼产生的升力为[19]
(7)
式中：Y为飞机升力(N)；
ρ为空气密度(kg/m3)；
S为机翼面积(m2)，通过查询飞机设计参数获得；
Cy为升力系数，通过查图表获得；
v为飞机滑行速度(m/s)。
在标准条件下空气密度约为1.29 kg/m3，当飞机滑行速度为20 m/s时，B737-800的升力计算值约为17 kN。
飞机-道面的耦合作用体现为：在道面平整度的影响下，飞机对道面产生冲击作用，使道面产生竖向位移，而竖向位移将进一步加强飞机振动，从而使飞机作用在道面上的力发生变化。通常跑道由面层、基层、土基组成，根据有限元方法，对道面结构离散后，可得道面的动力平衡方程见式(8)。
(8)
根据结构动力学振型叠加法[20]可知道面结构振动方程式(8)可表示成n阶振型的线性无关微分方程组。各阶模态均选取正交规格化的振型，可得到对应于第n阶振型的模态方程如式(9)。
(9)
式中：Pn1、Pn2为第n阶模态振型的广义力；
An、Bn为广义坐标；
ξn为第n阶模态的阻尼比；
ωn为第n阶模态的圆频率。
Pn1、Pn2包括飞机随机动荷载和静荷载，见式(10)
(10)
因飞机与道面始终接触，飞机-道面耦合状态下飞机轮胎竖向位移是道面振动产生的竖向位移与道面平整度之和，见式(11)。
(11)
式中：Zp5(x)、Zp6(x)为道面在左右轮胎位置处的平整度(m)；
Zr5(x)、Zr6(x)是道面模型中飞机左右轮胎作用位置处道面竖向振动位移(m)。
道面模型中道面任一横截面x处振动竖向位移可由各阶振型函数叠加而成，其表达式见式(12)。
(12)
式中：An、Bn与式(9)相同；
φn(x)为对应于道面某一水平位置x处的第n阶振型函数。
将式(12)代入式(11)，可得：
(13)
将飞机主起落架模型的振动平衡方程式(1)～(4)与道面动力平衡方程式(9)联立，并将式(13)代入，即可得到飞机-道面的耦合作用下的动力平衡方程组。方程组中有2N+4个方程式，N为道路子模型的模态方程组数，只取有代表性的几组低阶振型叠加即可。采用Newmark-β法可求解方程组，本文后续的研究中将通过有限元方法对动荷载求解。
由式(5)及式(6)可以看出，飞机滑行过程中主起落架上的动荷载主要受道面平整度、飞机升力的影响。以下探讨机场道面平整度的模拟方法。
采用滤波白噪声法可实现道面平整度的时域数值模拟。假设采用随机白噪声ω(t)作为系统输入，噪声单位强度为1，道面平整度位移函数q(t)为系统输出，则基于滤波白噪声的道面平整度一阶微分方程如式(14)表示[21]。由此看出，由滤波白噪声法模拟的道面平整度与道面参考空间频率功率谱密度Gq(n0)、飞机滑行速度v相关。
(14)
式中：nl为下截止空间频率，为0.08 m-1；
q(t)为道面平整度位移输入(m)；
ω(t)为单位强度为1的随机白噪声信号；
v为飞机滑行速度(m/s)；
n0为参考空间频率(m-1)；
Gq(n0)为参考空间频率n0的功率谱密度(m3)。
由于对道面功率谱密度(power spectral density, PSD)的测量相对较少，而国际平整度指数(international roughness Index，IRI)为我国机场必需通报的参数且测量相对便捷，因此可根据研究确定的IRI与PSD转换关系式(15)，计算得出不同IRI值对应的功率谱密度，见表1。
表1 IRI与PSD转换表Tab.1 Conversion between IRI and PSD
(15)
式中：a0为常数，a0=103m-1.5；
IRI为国际平整度指数(m/km)。
图2 道面平整度仿真程序Fig.2 Pavement flatness simulation program
当飞机滑行速度为20 m/s，道面平整度指数IRI分别为1.56 m/km、3.12 m/km、4.68 m/km时，利用图2中的仿真器确定道面平整度曲线见图3。由图3可看出，随着IRI的增大，平整度曲线的高差绝对值呈现增大趋势，且IRI不同时，同样长度的道面平整度变化趋势完全一致。
(a) IRI=1.56 m/km
(b) IRI=3.12 m/km
(c) IRI=4.68 m/km图3 不同IRI时道面平整度时域仿真曲线Fig.3 Time domain simulation curves of pavement flatness under different IRI
飞机动荷载的求解采用有限元方法，首先建立飞机四自由度有限元模型。由质量单元MASS21模拟飞机起落架悬挂部分、非悬挂部分的质量和转动惯量；
利用弹簧单元COMBIN14模拟悬挂部分与非悬挂部分之间的弹性连接(包括连接刚度和阻尼)；
利用连接单元MPC184模拟机身连接两侧主起落架的刚性梁。需同时为飞机悬挂系统质量单元赋予转动惯量属性，否则两个主起落架之间不能产生相互作用从而使有限元结果失真。
以常用民航客机B737-800为例，介绍由有限元法求解飞机动荷载的过程。此飞机含有一个前起落架，两个主起落架，其主起落架荷载分配系数为95%，前起落架仅承担5%的荷载重量，故采用四自由度主起落架模型。B737-800参数如表2所示[23]。
表2 B737-800飞机参数Tab.2 Parameters of B737-800
选择瞬态分析中完全法进行飞机随机荷载的求解。求解过程中，将道面平整度数值作为位移约束施加在机轮节点上，作用时间的长短由飞机的滑行速度决定，瞬态分析求解结束后即可得到道面不平整引起的起落架的随机动荷载数列。可得到不同机场道面平整度等级、不同滑行速度下飞机动荷载数列。当IRI为1.56 m/km，飞机以20 m/s匀速滑行时，B737-800的动荷载如图4所示。由图4可看出，受到升力的影响，飞机一个主起落架上的动荷载序列以280 kN(静荷载，飞机重量减去升力)为基准随机变化。
图4 飞机动荷载曲线Fig.4 Dynamic load curve of B737-800
(1) 不同IRI下飞机动荷载
IRI分别取3.12 m/km与4.68 m/km，假设飞机以20 m/s的速度匀速滑行，计算不同平整度等级下飞机每个主起落架上的动荷载，如图5所示。
图5 不同IRI下B737-800飞机动荷载Fig.5 Dynamic loads of B737-800 aircraft under different IRI
由图4及图5可以看出，当飞机滑跑速度不变，随着IRI增加，即跑道平整度等级由好变差时，飞机动荷载序列呈现整体增大的趋势。平整度指数为1.56 m/km、3.12 m/km、4.68 m/km时，B737-800每个飞机主起落架上的荷载峰值分别为290 kN、302 kN、311 kN，增长了7%。
(2) 不同滑行速度下飞机动荷载
飞机分别以10 m/s、20 m/s、30 m/s、40 m/s的速度匀速滑行于IRI为1.56 m/km的道面时，240 m的道面长度范围内，其动荷载变化曲线如图6所示。由图6可以看出，随着飞机滑跑速度的增大，飞机动荷载呈减小的趋势。究其原因为当飞机滑行速度超过10 m/s后，飞机升力对动荷载的影响大于平整度对动荷载的影响。飞机升力与滑跑速度呈正相关幂指数规律变化，当飞机滑跑速度在10～30 m/s时，飞机升力对飞机动荷载的作用逐渐显现，而飞机滑跑速度超过30 m/s后，飞机升力对飞机动荷载的影响明显加强，这与文献[24]研究的成果一致。
(a) v=10 m/s
(b) v=20 m/s
(c) v=30 m/s
(d) v=40 m/s图6 不同滑行速度下B737-800动荷载Fig.6 Dynamic load of B737-800 at different taxiing speeds
道面平整度激励下产生的飞机动荷载为随机函数，要更深入地研究飞机随机荷载，不仅要考虑其时域特性，更必须从频域着手对其变化规律进行探讨。本节采用Welch法分析飞机动荷载的功率谱密度，确定飞机动荷载的振动基频及其在频域内的能量分布，研究道面平整度、飞机滑行速度、悬挂部分的刚度及阻尼，非悬挂部分刚度及阻尼对飞机动荷载功率谱密度的影响规律。
Welch法又称为改进的平均周期图法，此方法采用信号重叠分段、加窗函数和FFT算法等计算一个信号序列的自功率谱估计[25]。Matlab中有固定的函数模块，可直接调用求解。
(1) 道面平整度对动荷载功率谱密度的影响
三种平整度条件下(IRI=1.56 m/km、IRI=3.12 m/km、IRI=4.68 m/km)飞机动荷载的功率谱密度如图7所示，功率谱密度峰值及对应的频率如表3所示。由图7与表3可知，飞机动荷载功率谱密度的峰值随道面等级的降低而显著增大，而峰值对应的频率不变。当IRI由1.56 m/km增大至4.68 m/km，即平整度状况由“好”变至“差”时，功率谱密度峰值由0.516 5×107N2·Hz-1增大至5.748×107N2·Hz-1，增大了10余倍。
图7 不同IRI时飞机荷载功率谱密度曲线Fig.7 Power spectral density curves of aircraft load under different IRI
表3 不同IRI时飞机动荷载基频Tab.3 Fundamental frequency of aircraft dynamic load with different IRI
(2) 飞机滑行速度对动荷载功率谱密度的影响
当IRI=1.56 m/km，飞机滑行速度为10 m/s、20 m/s、30 m/s与40 m/s时，飞机动荷载的功率谱密度如图8所示，功率谱密度峰值及峰值对应的频率如表4所示。分析图8与表4可知，随飞机滑跑速度增大，飞机动荷载功率谱密度频率亦增大，而频率对应的幅值先增大后减小。这说明随着滑跑速度的增加，动荷载的振动加强，对道面的破坏作用加大；
而幅值减小是受升力的影响，随着飞机滑跑速度逐渐增大，升力对动荷载的影响越来越显著，因此功率谱密度幅值呈现出先增大后减小的规律。
表4 不同滑行速度下的飞机动荷载基频Tab.4 Fundamental frequency of aircraft dynamic load at different taxi speeds
图8 不同滑行速度下飞机荷载功率谱密度曲线Fig.8 Power spectral density curves of aircraft load at different taxiing speeds
(3) 飞机悬挂系统参数对动荷载功率谱密度的影响
飞机以20 m/s的速度滑行于IRI为1.56 m/km的机场道面。调整飞机悬挂系统刚度，当悬架刚度由1.14×105N/m增大至6.14×105N/m，飞机动荷载功率谱密度峰值及峰值对应的频率如表5所示。分析表5可知：随着悬架刚度的增加，飞机动荷载功率谱密度的峰值增加，峰值的对应频率亦增大。这说明随着悬架刚度的增加，动荷载的振动加强，对道面的破坏作用加大。因此，进行飞机起落架设计时，在保证起落架承载能力及乘客乘坐舒适性前提下，应尽量选用小的起落架刚度，以减小飞机动荷载，提高道面的使用寿命。
表5 不同悬挂系统刚度时飞机动荷载基频Tab.5 Fundamental frequency of aircraft dynamic load with different suspension stiffness
其他参数不变，调整飞机悬挂系统的阻尼系数，阻尼由1×105N·s·m-1增加至7×105N·s·m-1，飞机动荷载功率谱密度峰值及峰值对应的频率如表6表示。分析表6可知：随着飞机悬挂部分阻尼的增加，动荷载功率谱密度的峰值降低，频率不变。由此可见，增加悬挂部分的阻尼可以减小飞机对跑道的动荷载，但悬挂系统阻尼对飞机动荷载能量分布的影响要小于悬挂系统的刚度。
表6 不同悬挂系统阻尼时飞机动荷载基频Tab.6 Fundamental frequency of aircraft dynamic load with different suspension system damping
(4) 飞机非悬挂系统参数对动荷载功率谱密度的影响
飞机以20 m/s的速度滑行于IRI值为1.56 m/km的机场道面。其他参数不变，调整飞机非悬挂系统的刚度，刚度由2×106N·m-1增加至6×106N·m-1，飞机动荷载功率谱密度峰值及峰值对应的频率如表7表示。分析表7可知：随着非悬挂系统刚度的增加，飞机动荷载功率谱密度的峰值增加，峰值对应的频率亦增大。这说明随着非悬架系统刚度的增加，动荷载的振动加强，对道面的破坏作用加大。而与悬挂系统刚度相比，非悬挂系统刚度对动荷载的影响偏小。
表7 不同非悬挂系统刚度时飞机动荷载基频Tab.7 Fundamental frequency of aircraft dynamic load with different stiffness of non-suspension system
其他参数不变，调整飞机非悬挂系统的阻尼系数，阻尼由2 000 N·s·m-1增加至10 000 N·s·m-1，飞机动荷载功率谱密度峰值及峰值对应的频率如表8所示。分析表8可知：随着飞机非悬挂部分阻尼的增加，动荷载功率谱密度的峰值稍有降低，频率不变。和其它参数相比，非悬挂部分的阻尼对动荷载的影响很小。
表8 不同非悬挂部分阻尼时飞机动荷载基频Tab.8 Fundamental frequency of aircraft dynamic load with different damping of non-suspension parts
建立了飞机四自由度模型，由滤波白噪声法模拟了道面平整度，在分析飞机-道面的耦合作用的基础上，将道面平整度作为约束施加至B737-800飞机四自由度模型，求解了不同平整度、不同滑行速度下飞机主起落架动荷载，并在时域与频域内分析了飞机随机动荷载的特征及变化规律。具体结论如下：
(1) 飞机滑跑速度不变，随着IRI增加，飞机随机动荷载呈现增大的趋势，平整度指数为1.56 m/km、3.12 m/km、4.68 m/km时，B737-800飞机每个主起落架上的动荷载峰值分别为290 kN、302 kN、311 kN，增长了7%，因此提高道面平度等级是减小飞机动荷载、提高道面使用寿命的有效措施之一。
(2) 当飞机滑跑速度由10 m/s增加至40 m/s，飞机动荷载减小，究其原因为机升力与滑跑速度呈幂函数规律变化，当飞机滑行速度超过10 m/s后，飞机升力对动荷载的影响远大于平整度对动荷载的影响。
(3) B737-800作用于道面上的随机动荷载以低频振动为主，能量主要集中在3 Hz以内；
随着飞机滑跑速度的提高，随机动荷载的振动频率增加，受升力的影响，动荷载功率谱密度峰值先升高后降低；
随着道面等级的降低，飞机随机动荷载功率谱密度的峰值增大，峰值对应的频率不变。
(4) 随着悬架与非悬挂系统刚度的增加，飞机随机动荷载功率谱密度的峰值增加，峰值的对应频率亦增大；
随着飞机悬挂部分阻尼的增加，随机动荷载功率谱密度的峰值明显降低，但频率保持不变；
而非悬挂部分阻尼，对动荷载功率谱密度的影响不明显。