1960年美国国有公路运输管理员协会(American association of state highway and transportation officials,AASHTO)道路试验研究表明，大约95％的路面服务性能来自于道路表面的平整度；长期路面性能(Long term pavement performance,LTPP)的研究也表明，路面平整度严重影响路面使用寿命^[1]。目前国际上应用最广的平整度评价指标是国际平整度指数(International roughness index,IRI)，美国和大部分欧洲国家均采用IRI为平整度验收指标^[2]。但该指标于机场道面平整度评价中的适用性问题越来越受到质疑^[2]，主要包括：不能反映飞机在不同轮迹带道面激励下飞机综合的振动反应；不能反映飞机所有机轮覆盖范围内的道面起伏情况对飞机滑行颠簸的影响效果；没有充分考虑飞机由于尺寸较大而与单条测线平整度相关性较差的情况。

基于此，本文将建立三自由度飞机力学模型，求解其振动方程，构造与该模型相应的平整度评价指标，并用Matlab软件里的Simulink工具对该模型进行建模并仿真，探究机场道面平整度的评价新方法。

1 模型建立及仿真

本文建立的三自由度飞机力学模型如图 1所示。图 1中：k₁₁，k₁₂，k₂₁，k₂₂，k₃₁和k₃₂分别为起落架前轮胎、前悬挂、左后轮胎、左后悬挂、右后轮胎和右后悬挂的刚度系数；c₁，c₂和c₃分别为前悬挂、左后悬挂和右后悬挂的阻尼系数；m₁，m₄，m₅和m₆分别为机身、前轮胎、左后轮胎和右后轮胎的质量；J_x,J_y分别为飞机模型绕x轴和y轴的转动惯量；q₁，q₂和q₃分别为前轮胎、左后轮胎和右后轮胎下的道面激励；Z₁，Z₄，Z₅和Z₆分别为机身、前轮胎、左后轮胎和右后轮胎的竖向位移；X，Y和Z分别为纵向(飞机滑行方向)、水平横向和竖向；θ和φ分别为飞机的横向滚转角和飞机的纵向俯仰角；O为飞机质心；a，b和2l分别为飞机质心到前起落架的距离、质心到两后起落架连线的垂直距离和两主起落架之间的距离。

本文建立的模型仿真系统包括3个部分，分别为道面激励子系统、振动处理子系统和指标处理子系统。

1.1 道面Simulink模型

大量实测数据表明道(路)面高程为平稳、遍历、均值为零的Gaussian过程，在时域内是各态历经的平稳随机过程^[3]。可以用不同形式的三角级数进行重构^[4]，因此将道面看做是一组波长、振幅和相位不同的正弦波的叠加，空间频率n(波长的倒数)范围为n₁,n₂，将区间n₁,n₂划分为N个小区间，第i个小区间的中心频率为n_mid-i，区间宽度为Δn，所对应的道面位移功率谱密度为G_q(n_i)，则

$\Delta n = {{{n_2} - {n_1}} \over N}$

(1)

${n_{mid - i}} = {n_1} + \left( {i - 1/2} \right)\Delta n$

(2)

${G_q}({n_i}) = {G_q}({n_{mid - i}}) = {G_q}({n_0}){({n_{mid - i}}{n_0})^{ - \omega }}$

(3)

式中：n₀为参考空间频率；ω为频率指数；G_q(n₀)为不平度系数，不同G_q(n₀)下的路面谱值见表 1^[5-6]；则该区间的功率谱为G_qn_iΔn，由巴什瓦等式^[7]得其标准差为G_qn_mid-iΔn，所以该区间对应的正弦信号为

$\sqrt {2{G_q}{n_{mid - i}}\Delta n} sin\left( {2\pi {n_{mid - i}}x + {\alpha _i}} \right)$

(4)

将N个这样的正弦波叠加起来得到道面随机激励为

${q_j}\left( x \right) = \sum\limits_{i = 1}^N {\sqrt {2{G_q}{n_{mid - i}}\Delta n} } sin\left( {2\pi {n_{mid - i}}x + {\alpha _i}} \right)$

(5)

式中：x为道面飞机滑行方向距起点的距离；q_j(x)为第j条轮迹带下面的道面激励(j=1，2，3) ；N为将频率范围n₁,n₂均分所得小区间数量；α_i为均匀分布于0,2π的随机数。

以文献^[8]中建立标准路面模型的方法为依据选择相关参数：n₀=0.1 m^-1；n₁=0.011 m^-1；n₂=2.830 m^-1；N=500；ω=2.00。

根据式(1~5) ，利用Matlab编写道面激励程序，可绘出3条A级道面激励沿道面的纵向分布，如图 2所示。

使道面等级分别为B，C，D，并根据表 1更换对应的G_q(n₀)的值，保持其他参数值不变，运行道面激励生成程序，可得相应等级道面激励的纵向分布，如图 3~5所示。

调用Simulink工具箱里面的相应模块建立三自由度飞机的机场道面平整度模型，模拟评价系统中的道面激励子系统，如图 6所示。该子系统包含模块：From Workspace1为从Matlab工作空间获取模块1；From Workspace2为从Matlab工作空间获取模块2；From Workspace3为从Matlab工作空间获取模块3；Gain 1为增益模块1；Gain 2为增益模块2；Gain 3为增益模块3；Add为加模块；Out为输出端口模块。其中：From space模块的功能为从工作空间和存储工作区中读取数据作为输入信号；Gain模块的功能为使输入的信号乘以一个向量，使信号转化为向量的形式；Add模块的作用是对两个或多个信号进行求和运算；Out模块作用是将信号输入到下一级。

1.2 三自由度飞机的动力学、数值及Simulink模

建立飞机三自由度物理模型如图 1，模型中机体为刚体，不考虑自身变形对其运动状态的影响。

由于本文研究的道面类型主要为水泥混凝土道面，刚度很大，变形很小，可以不用考虑飞机轮胎和道面的刚度耦合作用^[9]，且轮胎刚度系数对振动的影响要远大于其阻尼系数的影响^[10]，因而本文模型不计轮胎阻尼，对于图 1所示的动力学模型，利用拉格朗日方程建立其振动微分方程为

$\left[ M \right]\{ \ddot Z\} + \left[ C \right]\{ \dot Z\} + \left[ K \right]\left\{ Z \right\} = [{K_t}]\left\{ q \right\}$

(6)

式中

$\left\{ Z \right\} = {\{ {z_1},\varphi ,\theta ,{z_4},{z_5},{z_6}\} ^T}$

(7)

$\{ \dot Z\} = {\{ {{\dot z}_1},\dot \varphi ,\dot \theta ,{{\dot z}_4},{{\dot z}_5},{{\dot z}_6}\} ^T}$

(8)

$\{ \ddot Z\} = {\{ {{\ddot z}_1},\ddot \varphi ,\ddot \theta ,{{\ddot z}_4},{{\ddot z}_5},{{\ddot z}_6}\} ^T}$

(9)

$\left\{ q \right\} = {\{ {q_1},{q_2},{q_3}\} ^T}$

(10)

${\left[ M \right]_{6 \times 6}} = diag\{ {m_1},{J_y},{J_x},{m_4},{m_5},{m_6}\} $

(11)

${[{K_t}]_{6 \times 3}} = \left\lfloor {\matrix{ 0 & 0 & 0 \cr 0 & 0 & 0 \cr 0 & 0 & 0 \cr {{k_{11}}} & 0 & 0 \cr 0 & {{k_{21}}} & 0 \cr 0 & 0 & {{k_{31}}} \cr } } \right\rfloor $

(12)

${\left[ C \right]_{6 \times 6}} = \left\lfloor {\matrix{ {{c_1} + {c_2} + {c_3}} & { - {c_1}a + {c_2}b + {c_3}b} & 0 & { - {c_1}} & { - {c_2}} & { - {c_3}} \cr { - {c_1}a + {c_2}b + {c_3}b} & {{c_1}{a^2} + {c_2}{b^2} + {c_3}{b^2}} & 0 & {{c_1}a} & { - {c_2}b} & { - {c_3}b} \cr 0 & 0 & {{c_2}{l^2} + {c_3}{l^2}} & 0 & {{c_2}l} & { - {c_3}l} \cr { - {c_1}} & {{c_1}a} & {{c_1}} & 0 & 0 & 0 \cr { - {c_2}} & { - {c_2}b} & {{c_2}l} & 0 & {{c_2}} & 0 \cr { - {c_3}} & { - {c_3}b} & { - {c_3}l} & 0 & 0 & {{c_3}} \cr } } \right\rfloor $

(13)

${\left[ K \right]_{6 \times 6}} = \left\lfloor {\matrix{ {{k_{12}} + {k_{22}} + {k_{32}}} & { - {k_{12}}a + {k_{22}}b + {k_{32}}b} & 0 & { - {k_{12}}} & { - {k_{22}}} & { - {k_{32}}} \cr { - {k_{12}}a + {k_{22}}b + {k_{32}}b} & {{k_{12}}{a^2} + {k_{22}}{b^2} + {k_{32}}{b^2}} & 0 & {{k_{12}}a} & { - {k_{22}}b} & { - {k_{32}}b} \cr 0 & 0 & {{k_{22}}{l^2} + {k_{32}}{l^2}} & 0 & {{k_{22}}l} & { - {k_{32}}l} \cr { - {k_{12}}} & {{k_{12}}a} & 0 & {{k_{12}} + {k_{11}}} & 0 & 0 \cr { - {k_{22}}} & { - {k_{22}}b} & {{k_{22}}l} & 0 & {{k_{22}} + {k_{21}}} & 0 \cr { - {k_{32}}} & { - {k_{32}}b} & { - {k_{32}}l} & 0 & 0 & {{k_{32}} + {k_{31}}} \cr } } \right\rfloor $

(14)

根据式(6~14) ，利用Simulink工具箱中Continuous模块库中的相关模块建立模型系统的振动处理子系统，如图 7所示。该系统所包含模块：In为输入模块；Gain 1为增益模块1；Gain 2为增益模块2；Gain 3为增益模块3；Gain 4为增益模块4；Gain 5为增益模块5；Gain 6为增益模块6；Gain 7为增益模块7；Add为加模块；Integrator 1为积分模块1；Integrator 2为积分模块2；Out为输出端口模块。其中：In模块的作用是做为子系统里面的一个接受外部输入的一个端口从子系统的上一级接受变量向量；Integrator模块的作用是对信号进行积分。

三自由度飞机模型可选取任意一种起落架形式为前三点式飞机的参数做为计算参数，本文以波音707-320C型飞机原始数据^[11-13]为例，选取模型参数如下：m₁=1.44×10⁵ kg;m₄=156 kg;m₅=m₆=2294 kg;J_x=5.28×10⁶ kg·m²；J_y=7.28×10⁶ kg·m²;a=16.66 m;b=1.32 m;2l=6.73 m;k₁₁=2.37×10⁶ N/m;k₂₁=k₃₁=1.7×10⁷ N/m;k₁₂=1.92×10⁶ N/m;k₂₂=k₃₂=5.54×10⁶ N/m;c₁=3.61×10⁴ (N·(m/s)^-1);c₂=c₃=2.13×10⁵ (N·(m/s)^-1)。

将以上参数代入式(11~14) 得

${\left[ M \right]_{6 \times 6}} = \left\lfloor {\matrix{ {1.44 \times {{10}^5}} & 0 & 0 & 0 & 0 & 0 \cr 0 & {7.28 \times {{10}^6}} & 0 & 0 & 0 & 0 \cr 0 & 0 & {5.28 \times {{10}^6}} & 0 & 0 & 0 \cr 0 & 0 & 0 & {156} & 0 & 0 \cr 0 & 0 & 0 & 0 & {2294} & 0 \cr 0 & 0 & 0 & 0 & 0 & {2294} \cr } } \right\rfloor $

(15)

${\left[ C \right]_{6 \times 6}} = \left\lfloor {\matrix{ {4.62 \times {{10}^5}} & { - 3.91 \times {{10}^4}} & 0 & { - 3.61 \times {{10}^4}} & { - 2.13 \times {{10}^5}} & { - 2.13 \times {{10}^5}} \cr { - 3.91 \times {{10}^4}} & {1.08 \times {{10}^7}} & 0 & {6.01 \times {{10}^5}} & { - 2.81 \times {{10}^5}} & { - 2.81 \times {{10}^5}} \cr 0 & 0 & {4.82 \times {{10}^6}} & 0 & {7.17 \times {{10}^5}} & { - 7.17 \times {{10}^5}} \cr { - 3.61 \times {{10}^4}} & {6.01 \times {{10}^5}} & 0 & {3.61 \times {{10}^4}} & 0 & 0 \cr { - 2.13 \times {{10}^5}} & { - 2.81 \times {{10}^5}} & {7.17 \times {{10}^5}} & 0 & {2.13 \times {{10}^5}} & 0 \cr { - 2.13 \times {{10}^5}} & { - 2.81 \times {{10}^5}} & { - 7.17 \times {{10}^5}} & 0 & 0 & {2.13 \times {{10}^5}} \cr } } \right\rfloor $

(16)

${{\left[ K \right]}_{6\times 6}}=\left\lfloor \begin{matrix} 1.3\times {{10}^{7}} & -1.74\times {{10}^{7}} & 0 & -1.92\times {{10}^{6}} & -5.54\times {{10}^{6}} & -5.54\times {{10}^{6}} \\ -1.74\times {{10}^{7}} & 5.52\times {{10}^{8}} & 0 & 3.20\times {{10}^{7}} & -7.31\times {{10}^{6}} & -7.31\times {{10}^{6}} \\ 0 & 0 & 1.25\times {{10}^{8}} & 0 & 1.86\times {{10}^{7}} & -1.86\times {{10}^{7}} \\ -1.92\times {{10}^{6}} & 3.20\times {{10}^{7}} & 0 & 4.29\times {{10}^{6}} & 0 & 0 \\ -5.54\times {{10}^{6}} & -7.31\times {{10}^{6}} & 1.86\times {{10}^{7}} & 0 & 2.25\times {{10}^{7}} & 0 \\ -5.54\times {{10}^{6}} & -7.31\times {{10}^{6}} & -1.86\times {{10}^{7}} & 0 & 0 & 2.25\times {{10}^{7}} \\ \end{matrix} \right\rfloor $

(17)

${\left[ {{K_t}} \right]_{6 \times 3}} = \left\lfloor {\matrix{ 0 & 0 & 0 \cr 0 & 0 & 0 \cr 0 & 0 & 0 \cr {2.37 \times {{10}^6}} & 0 & 0 \cr 0 & {1.7 \times {{10}^7}} & 0 \cr 0 & 0 & {1.7 \times {{10}^7}} \cr } } \right\rfloor $

(18)

将式(15~18) 中的数值输入到图 7所示的振动处理子系统中，设置道面等级依次为A，B，C和D，运行系统，得到与图 3~6对应的飞机振动响应量Z₁，aφ和lθ，如图 8~11所示。

1.3 新评价指标FARI的构建

由于IRI只能反映1/4车模型悬挂的竖向累积位移，不能反映飞机滑行过程中对机械磨损、乘客舒适性和燃油消耗等影响较为显著的飞机横向的滚转效应和纵向的俯仰效应。

IRI表达式^[14]为

$IRI=\frac{1}{L\prime }_{0}^{L\prime }\left| Z_{2}^{'}-Z_{1}^{'} \right|dx$

(19)

IRI来源的1/4车模型如图 12所示。

图中：L^'为IRI输出间隔长度；Z^'₀ 为道面不平度；Z^'₁ 为1/4车模型轮胎竖向位移；Z^'₂ 为1/4车模型车体坚向位移；k^'₁ 为1/4车模型轮胎刚度系数；k^'₂ 为1/4车模型悬架刚度系数；c'为1/4车模型悬架阻尼系数；m^'₁ 为1/4车模型非悬挂系统质量(车轮质量)；m^'₂ 为1/4车模型悬挂系统质量(车身质量)。

因此，依据1/4车模型推导得IRI的原理，由三自由度飞机模型推得整机平整度指数(Full aircraft roughness index,FARI)，表达式为

$FARI = {1 \over L}\left( {_0^L\left| {{Z_1}} \right|dx + _0^L\left| {a\varphi } \right|dx + _0^L\left| {l\theta } \right|dx} \right)$

(20)

式中：L为FARI输出间隔长度；Z₁为飞机模型质心处纵向位移；a为飞机模型质心到前起落架的距离；l为飞机模型主起落架间距的一半；φ为飞机模型俯仰角；θ为飞机模型滚转角。

由式(20) 可知，新指标FARI由3部分组成，分别为：飞机质心处的累积竖向位移量、俯仰效应引起的前起落架处机身累积位移量以及滚转效应引起的位于飞机质心处机身横截面最边缘位置的累积竖向位移量。

根据式(20) ，Simulink工具箱Continuous模块库中的相关模块建立模型系统的指标处理子系统如图 13所示。该子系统中所包含模块：In为信号输入模块；Abs为取绝对值模块；Gain为增益模块；Pulse generator为脉冲发生器模块；Integrator limited为定积分模块；Out为输出端口模块。其中：Abs模块的作用是对信号取绝对值使其全部为正；Pulse generator模块的作用是以一定的间隔产生脉冲信号；Integrator limited模块的作用是以脉冲信号间隔为界限对信号进行定积分。

道面等级依次取A，B，C和D，运行系统，得到不同道面等级的FARI值沿道面纵向分布如图 14~17所示。

由图 14~17可知：A级道面激励下整机平整度指数取值范围为0~4 m/km；B级道面激励下整机平整度指数取值范围为0~8 m/km；C级道面激励下整机平整度指数取值范围为0~20 m/km；D级道面激励下整机平整度指数取值范围为0~46 m/km。由图知FARI的值随道面等级的降低而升高，这是由于随着道面等级的降低道面越来越不平整，从而引起的飞机颠簸量越来越大，单位距离内的飞机累积颠簸量也越来越大，所以FARI值越来越大，这与实际情况相符合。

2 评价指标分析

目前平整度评价指标较多，对平整度评价指标进行评判的标准一般从以下4个方面考虑，即飞机的油耗、飞机的机械磨损、道面的损坏程度或者乘客的舒适性。本文将从乘客的舒适性方面考虑，通过研究FARI与其相关性，从而对FARI进行评判。目前国际上使用最广泛的振动舒适性标准是ISO 2631-1—1997^[15]和BS 6841—1987^[16]，两个标准均采用加权加速度均方根值作为基本评价指标^[2]，其计算公式如下

${a_w} = {\left[ {{1 \over X}_0^Xa_W^2\left( x \right)dx} \right]^{{1 \over 2}}}$

(21)

式中：X为分析距离长度；a_w(x)为瞬时频率加权加速度幅值。

通过研究FARI，IRI和a_w(x)的相关性，得到更为优越的道面平整度评价指标。相关系数计算公式如下

$r = {{n\sum\limits_{i = 1}^n {{u_i}{v_i}} - \sum\limits_{i = 1}^n {{u_i}} \times \sum\limits_{i = 1}^n {{v_i}} } \over {\sqrt {n\sum\limits_{i = 1}^n {u_i^2 - {{\left( {\sum\limits_{i = 1}^n {{u_i}} } \right)}^2}} } \times \sqrt {n\sum\limits_{i = 1}^n {v_i^2} - {{\left( {\sum\limits_{i = 1}^n {{v_i}} } \right)}^2}} }}$

(22)

式中：n为所取数据的组数；u_i为整机平整度指数FARI或3条测线的国际平整度指数IRI的平均值；v_i为加速度均方根值。

将各子系统系统进行组合得完整仿真系统如图 18所示。系统中所包含模块：Subsystem 1为道面激励子系统模块；Subsystem 2为振动处理子系统模块；Subsystem 3为指标处理子系统模块；Scope为示波器模块；Out为输出端口模块。

模拟A~D'八个等级的8条跑道道面的高程数据，其中A'，B'，C'和D'为介于上下相邻两个道面等级之间的道面等级，其值为上下相邻两个道面等级不平度系数值的中值。模拟道面长度4 000 m，离散间隔1 m，飞机参数参考美国波音707-320C型飞机，运行基于三自由度飞机模型的道面平整度评价仿真系统，得到相应的整机平整度指数FARI，通过Matlab编程，计算对应不同等级道面的飞机的加权加速度均方根值a_w(该加速度为飞机质心处竖向加速度、俯仰效应引起的前起落架处机身竖向加速度机身以及滚转效应引起的位于飞机质心处机身横截面最边缘位置的竖向加速度三者的均值)，结果见表 2。

依据表 2数据，由式(21) 计算的整机平整度指数FARI和加权加速度均方根值a_w的相关系数r₁=0.980 6。将a_w作为因变量y，FARI作为自变量x，得拟合表达式为

$y = 0.043 8x + 0.484 3$

(23)

参照世界银行公布的IRI计算源程序，用Matlab编程计算A~D'八个等级道面的IRI值，取每个等级道面3条测线的IRI平均值为该级道面的IRI值，并计算对应不同等级道面的飞机的加权加速度均方根值a_w，结果见表 3。

依据表 3数据，由式(22) 计算的国际平整度指数IRI和加权加速度均方根值a_w的相关系数r₂=0.886 9。对比可知r₁＞r₂，即FARI与a_w的相关程度要大于IRI与a_w的相关程度，因此用FARI来表示道面的平整度要比IRI更为合理。

3 评价指标的量化分级

由于整机平整度指数FARI与加权加速度均方根值a_w有着较好的相关性，因此可以依据a_w的标准对FARI进行量化分级。加权加速度均方根值与行车舒适性的关系见表 4^[17]。

对比表 2，4，以表 4中人的舒适度标准为依据，可以将道面不平度分为3个等级：当a_w=0.63时，由式(23) 整机平整度指数FARI=3.3，此时的加速度值是舒适与不太舒适的界限；当a_w=1.25时，得整机平整度指数FARI=17.5，此时的加速度值是不舒适与很不舒适的界限。以这两个界限为分界点对FARI分级如表 5所示。

表 5即为整机平整度指数的分级标准，对于波音707-320C型飞机，当FARI值小于3.3时，道面平整度的情况为良好，当FARI大于等于3.3且小于等于17.5时，道面平整度的情况为中等，当FARI大于17.5时，道面平整度的情况为差。

4 结束语

本文构建了道面激励数学模型，建立了飞机的三自由度力学振动模型，给出并求解了其振动方程，提出了新的平整度评价指标FARI。通过研究发现FARI比IRI更适合于机场道面平整度评价，并对FARI进行了分级，给出了分级标准。本文结果为机场道面平整度评价提供了新视角和新方法。