国内外均开展了针对结冰后飞机的飞行动力学特性的研究，主要集中在对飞机的总体性能、操纵性、稳定性、自动驾驶仪工作状态下的闭环特性等方面。NASA早于20纪80年代开展了一系列的飞机结冰研究计划，Cebeci[2]、Bragg[3-7]、Khodadoust[8-10]、Ranaudo[11-12]和Potapczuk[13-16]等取得了一些重要的研究成果。Hiltner[17]基于飞行试验数据与风洞试验数据，系统地研究了平尾结冰情形下的动力学特性，指出平尾结冰会使得飞机出现明显的稳定性与操纵性变差的问题。Sibilski等[18]对飞机爬升时遭遇结冰情形下的动力学响应进行了分析，得出当飞机从水平飞行转为爬升时，由于冰形的积聚容易引发飞机失速，并且指出在失速点附近，由于飞机的法向过载接近于1，机组人员很难意识到飞机会突然失控的危险。Cunningham[19]在建立的简化的结冰模型基础上，较为系统地分析了飞机遭遇结冰后，飞机操纵性、配平特性、纵横向模态特性等动力学性能的改变。中国对于结冰后飞行动力学特性的研究起步相对较晚，针对特定型号的飞机如Y12、Y7-200A和ARJ-21-700开展过结冰飞行试验，研究了真实环境下结冰对飞行性能和操稳特性的影响。北京航空航天大学[20-21]、南京航空航天大学[22-24]、中国空气动力研究与发展中心[25]等科研院所均开展了结冰对飞行动力学特性的研究。
随着中国民用大飞机和军用大型运输机的快速发展，基于现有有限认知能力发展起来的结冰防护技术，还不能完全消除结冰对飞行安全的危害，结冰条件下的安全保障需求迫切。为了较为全面地了解结冰对大飞机操稳特性的影响，本文基于数值模拟手段获得了背景飞机在干净外形、中度结冰和重度结冰3种情况下的气动数据库，构建了飞机六自由度非线性动力学模型。在仿真模型的基础上设计了俯仰角保持(Pitch Attitude Hold，PAH)、滚转角保持(Roll Attitude Hold，RAH)、高度保持(Altitude Hold，ALH)几种模式的自动驾驶仪。与国内外研究相比，重点针对结冰条件下飞机在开环和闭环下的操稳特性进行分析。通过开环仿真，研究积冰对飞机配平特性、纵向长短周期模态、横航向模态的影响，模拟飞机在积冰过程中的飞行状态和响应特性的变化。通过闭环仿真，研究积冰对自动驾驶仪工作性能的影响。本文研究成果可为后续开展的研究工作，如结冰情形下的安全操纵策略、结冰飞行风险评估及飞行安全边界保护等提供理论支撑。
1 背景飞机结冰气动数据的获取 目前，常用的飞机结冰气动模型主要采用伊利诺伊州大学Bragg教授提出的一种结冰参量模型[26]，该模型采用结冰因子来评估结冰对气动性能的影响，该参数与飞机自身的尺寸、飞行状态或飞机受结冰影响的敏感性相关，但该模型无法直接移植到其他气动布局差异较大的飞机上来估算结冰的气动影响。
为了研究典型大飞机机翼前缘积冰对飞机飞行动力学特性的影响，必须获得可靠的结冰气动数据。本文采用高精度的数值模拟方法获取飞机结冰气动参数。首先，构建了背景飞机的几何模型，建立了干净构型的三维数模，如图 1所示。考虑到角状冰对飞机气动特性影响最大，根据可靠的结冰冰形的试验数据形成了典型结冰环境下的机翼前缘重度和中度角状冰模型，如图 2所示。三维冰形从机翼翼根至翼稍沿弦长方向剖面形状如图 3所示。然后，基于高精度的数值模拟方法计算了典型状态下机翼结冰对飞机气动参数的影响。最终，通过风洞试验虚拟飞行验证，获得了一套可靠的飞机结冰后的气动参数数据库。
 |
| 图 1 背景飞机三维模型 Fig. 1 3D model of background aircraft |
| 图选项 |
 |
| 图 2 重度结冰和中度结冰几何模型对比 Fig. 2 Comparison of heavy and moderate icing geometric models |
| 图选项 |
 |
| 图 3 机翼前缘重度结冰的冰形剖面 Fig. 3 Ice shape profile of heavy icing on leading edge of wing |
| 图选项 |
图 4~图 6为飞机在3种不同结冰构型下(干净外形构型clean、中度结冰构型ice-B、重度结冰构型ice-A)，马赫数Ma=0.4，平尾安装角为0°，升降舵偏角为0°时的升力系数CL、阻力系数CD与俯仰力矩系数Cm随迎角α的变化曲线。图 4~图 6结果显示，干净外形、中度结冰和重度结冰状态下，升力线斜率和最大升力系数依次下降，失速迎角依次减小，最小阻力依次增加。在小迎角线性区内，纵向静稳定裕度依次降低，在大迎角区域呈现非线性的变化趋势。
 |
| 图 4 升力系数变化曲线 Fig. 4 Variation curves of lift coefficient |
| 图选项 |
 |
| 图 5 阻力系数变化曲线 Fig. 5 Variation curves of drag coefficient |
| 图选项 |
 |
| 图 6 俯仰力矩系数变化曲线 Fig. 6 Variation curves of pitching moment coefficient |
| 图选项 |
2 飞机的动力学模型 飞机本体六自由度全量的非线性动力学模型可以用式(1)的向量形式进行表示[27-28]：
 | (1) |
式中：u为控制向量；x为状态向量。
 | (2) |
式中：V、α和β分别为飞行速度、飞机迎角和侧滑角；p、q和r分别为滚转、俯仰和偏航角速度；q0、q1、q2、q3为四元素；xg、yg、zg为飞机在地面坐标系下的空间位置。
 | (3) |
式中：δth为飞机油门偏度；δe、δa和δr分别为升降舵、副翼和方向舵舵面的偏角。
3 自动飞行控制律模型 本文在积冰飞机动力学模型的基础上，设计具有俯仰角保持、滚转角保持、高度保持3种常用模式的自动飞行控制系统。
俯仰角保持模式主要功能为保持飞机给定的俯仰姿态角θref，其工作原理如图 7所示。
 |
| 图 7 俯仰角保持模式结构框图 Fig. 7 Structure diagram of pitch attitude hold mode |
| 图选项 |
通过增加俯仰角速率反馈作为控制增稳，增加俯仰角保持下飞机的稳定性和阻尼。通过俯仰角指令与飞机当前俯仰角反馈的差值Δθ，控制飞机跟踪保持俯仰姿态指令。控制律表示为
 | (4) |
式中：kθP、kθI、kq为控制增益。
滚转角保持模式主要功能为保持飞机给定的滚转姿态角?ref。在滚转角保持模式下，飞机通过对副翼的控制使飞机跟踪滚转角期望值。由于飞机的横向和侧向存在耦合，飞机滚转角保持要通过横向和航向2个通道进行调节。飞机侧向转弯一般采用无侧滑的协调转弯方式，满足：
 | (5) |
式中：?为滚转角; ψ为偏航角；g为重力加速度。由于滚转角?的存在，还要保证飞机垂直方向上的升力分量与重力平衡，水平方向上的分量与飞机的离心力平衡。工作原理如图 8所示。
 |
| 图 8 无侧滑滚转角保持模式结构框图 Fig. 8 Structure diagram of roll attitude hold mode without sideslip |
| 图选项 |
其控制律可以表示为
 | (6) |
式中：Δ?为飞机滚转角指令与飞机当前滚转角的差值，即Δ?=?ref－? ；k?P、k?I、kp、kψP、kr、kβ为控制增益。
高度保持模式是飞机控制系统中重要的组成部分，如图 9所示。图中飞机纵向内回路即为俯仰角保持回路。在飞机的爬升、巡航和着陆初始阶段都需要保持高度的稳定。高度保持模式下，飞机能自动维持在某一固定高度。通过直接引入期望高度Href，得到高度差信号，经过俯仰姿态保持内回路，控制飞机姿态变化，实现对飞行高度的控制。其控制律可以表示为
 | (7) |
 |
| 图 9 高度保持模式结构框图 Fig. 9 Structure diagram of altitude hold mode |
| 图选项 |
式中：ΔH为飞机高度指令与飞机当前高度的高度差，即ΔH=Href－H；kHP、kHI、kHD为控制增益。
4 仿真分析 4.1 本体操稳特性分析
4.1.1 不同程度积冰对飞机配平特性的影响 设定初始飞行状态为：高度3000m，速度120m/s，定直平飞。分别在干净外形、机翼前缘中度结冰和重度结冰3种情况下对飞机的配平参数进行比较，如表 1所示。
表 1 不同积冰情形下的飞机配平参数 Table 1 Aircraft trim parameters under different ice accretion conditions
| 积冰情形 | 配平迎角/(°) | 升降舵偏角/(°) | 油门偏度/% |
| 干净外形 | 6.83 | -4.37 | 29.36 |
| 中度结冰 | 8.03 | -2.17 | 36.90 |
| 重度结冰 | 8.67 | -1.84 | 45.73 |
表选项
对比发现，在相同的飞行状态配平条件下，由于飞机结冰严重程度不断增大，导致飞机的升力系数不断下降，阻力系数不断增大，因此所需要的配平迎角不断增大，中度结冰与重度结冰情况下平飞需要的发动机推力较干净外形增大了25.7%和55.8%。同时在小迎角线性区，机翼前缘积冰导致俯仰配平力矩减小，所需配平升降舵偏角绝对值随结冰严重程度的增大而减小。
4.1.2 不同程度积冰对飞机本体纵向长短周期模态的影响 在相同的飞行状态下，对不同积冰情形飞机本体纵向长短周期模态特征参数进行分析，如表 2所示。表中：ωnsp、ωnp分别为短、长周期模态自由振荡频率，ζsp、ζp分别为短、长周期模态的阻尼比，?nz/?α为加速度敏感性，CAP为操纵期望参数。
表 2 不同积冰情形下的飞机本体纵向长短周期模态特征参数 Table 2 Longitudinal long- and short-period modal characteristic parameters of aircraft body under different ice accretion conditions
| 积冰情形 | ωnsp/(rad·s-1) | ζsp | ωnp/(rad·s-1) | ζp |  | CAP/(rad·s-2·g-1) |
| 干净外形 | 1.7775 | 0.4566 | 0.0822 | 0.0282 | 7.3921 | 0.4274 |
| 中度结冰 | 1.4944 | 0.4962 | 0.1031 | 0.0381 | 6.0303 | 0.3703 |
| 重度结冰 | 1.5219 | 0.4907 | 0.0957 | 0.0557 | 6.4009 | 0.3618 |
表选项
可以看到，随着积冰严重程度的增大，短周期阻尼比ζsp先增大后减小，总体影响不大，但总体来说，机翼结冰相较于干净外形降低了纵向短周期模态自由振荡频率ωnsp。对于长周期模态来说，阻尼比是逐渐增大的，但影响也不大。操纵期望参数CAP则随着积冰严重程度的增大而逐渐减小，但变化幅度不大，根据MIL-F-8785C[29]可知，巡航阶段下这3种积冰情形下的CAP均处于1级飞行品质。
4.1.3 不同程度积冰对飞机本体横航向模态的影响 在相同的飞行状态下，对不同积冰情形飞机本体横航向模态的特征参数进行分析，如表 3所示。表中：τS为螺旋模态的时间常数。
表 3 不同积冰情形下的飞机本体横航向模态特征参数 Table 3 Lateral modal characteristic parameters of aircraft body under different ice accretion conditions
| 积冰情形 | τR/s | λS | τS/s | ωnDR/(rad·s-1) | ζDR | τDR/s |
| 干净外形 | 0.1551 | -0.0018 | 553.7114 | 1.2671 | 0.1965 | 4.0157 |
| 中度结冰 | 0.1636 | 0.0066 | 150.8256 | 1.3061 | 0.2137 | 3.5833 |
| 重度结冰 | 0.1676 | 0.0019 | 540.3449 | 1.2665 | 0.2123 | 3.7196 |
表选项
可以看到，对于滚转模态，时间常数τR随着积冰严重程度的增大而逐渐增大，说明积冰导致滚转模态收敛时间增大。对于螺旋模态，积冰后对应螺旋模态的特征值λS均大于0，说明螺旋模态缓慢发散，同时随着积冰严重程度的增大而加剧。对于荷兰滚模态，积冰导致荷兰滚模态阻尼比ζDR增大，时间常数τDR减小，自由振荡频率ωnDR中度结冰时变大，重度结冰时反而变小，说明影响程度与结冰严重程度没有明显的正相关性。
4.1.4 不同程度积冰对飞机响应特性的影响 为研究飞机不同程度结冰后的响应，设置初始配平的飞行状态不变，在t=10s时升降舵施以持续2 s的脉冲激励，升降舵偏角增量为-2°。通过动力学仿真，分别计算出飞机在干净外形、中度结冰与重度结冰3种情形下的动态响应，仿真结果如图 10所示。
 |
| 图 10 主要飞行参数响应曲线 Fig. 10 Response curves of main flight parameters |
| 图选项 |
从图 10中可以看到，对于长周期模态参数(高度H、速度V)而言，相较于干净外形情况，结冰后这些参数在波动过程中的峰值有所减小；对于短周期模态参数(迎角α、俯仰角速度q)来说，3种情况下的动态响应特征基本一致，但在积冰情形下由于失速临界迎角的降低，飞机在该升降舵激励下，响应过程中的最大迎角可能已处在失速迎角范围内，会使飞机的安全性受到威胁。
4.2 闭环系统仿真 为了分析遭遇中度结冰和重度结冰情形下对飞机自动驾驶仪性能的影响，分别在干净外形、中度结冰和重度结冰3种情况下进行仿真，对比仿真结果，分析积冰对自动驾驶仪工作性能的影响及飞机动态响应的变化。仿真初始条件均设定为3种情形下的配平量，总仿真时间设定为100 s。
4.2.1 俯仰角保持模式 在俯仰角保持模式下，t=0 s时令俯仰角保持在10°，仿真结果如图 11所示。
 |
| 图 11 俯仰角保持模式仿真结果 Fig. 11 Simulation results in pitch attitude hold mode |
| 图选项 |
通过仿真结果可以看到，干净外形和积冰情况下俯仰角跟随效果均比较理想，说明设计的俯仰角保持控制律的正确性。积冰情况下飞机迎角与干净外形情形相比变化不大，速度下降幅度较干净外形小。同时积冰情况下高度增长更为平缓，并随着结冰严重程度的增大，增长幅度越小。总体来看，在该状态下积冰的影响较小。
4.2.2 滚转角保持模式 在滚转角保持模式下，t=0 s时令滚转角保持在20°，仿真结果如图 12所示。
 |
| 图 12 滚转角保持模式仿真结果 Fig. 12 Simulation results in roll attitude hold mode |
| 图选项 |
从图 12中可以看出，干净外形和积冰情况下滚转角跟随效果均比较理想，说明设计的滚转角保持控制律的正确性。在滚转角保持模式下，相较于干净外形情况，积冰情况下飞行速度的振荡幅度较大，且重度结冰最严重，高度下降幅度最大。总体来看，在该状态下积冰的影响较小。
4.2.3 高度保持模式 在高度保持模式下，t=0s时令高度保持在3100m，仿真结果如图 13所示。
 |
| 图 13 高度保持模式仿真结果 Fig. 13 Simulation results in attitude hold mode |
| 图选项 |
从图 13可以看出，干净外形情况下高度的跟随效果均比较理想，说明设计的高度保持自动驾驶仪模型的正确性。相较于干净外形情况，积冰导致高度响应超调量增大。迎角响应在积冰与干净外形情形下最终都趋于平稳，积冰情形下迎角稳态值大于干净外形情形。速度在干净外形时逐渐趋于平稳，而在积冰情形下速度值缓慢下降，重度结冰时下降幅度最大。总体来看，在该状态下积冰的影响较小。
5 结论 本文着重研究大飞机机翼前缘积冰对操稳特性的影响。根据实验数据构建了典型的机翼前缘积冰冰形，通过数值模拟得到背景飞机的结冰气动数据，建立了飞机积冰动力学模型，设计了俯仰角保持、滚转角保持及高度保持的自动驾驶仪模型，通过开环与闭环仿真研究了积冰对飞机操纵的影响，在本文算例条件下得到以下结论：
1) 积冰易导致飞机配平特性的恶化，主要体现在：积冰下配平需要更大的迎角和发动机推力，这是由于机翼结冰减小了飞机升力系数，增大了阻力系数。
2) 积冰将降低纵向短周期模态自由振荡频率，增大长短周期的阻尼比，并减小操纵期望参数CAP，但总体影响较小。同时积冰将导致滚转模态收敛时间增大；积冰会使螺旋模态发散，荷兰滚模态变化不大。
3) 积冰对自动驾驶飞机指令的跟踪总体影响较小，最主要的变化体现在高度与速度上。
4) 针对本文中的算例，对操稳特性的影响程度与机翼结冰严重程度没有明显的正相关性，这是由于结冰增长带来的气动力变化的非线性导致的。
本文主要研究线性区迎角范围内的积冰对飞机操稳特性的影响。下一步将研究积冰情形下大迎角飞行动力学特性及积冰增长导致的非定常空气动力学与飞行力学的非线性耦合规律，为研究积冰动态增长过程中的空气动力学特性和飞行力学特性奠定基础。
参考文献
| [1] | REEHORST A L, JR ADDY H E, COLANTONIO R O.Examination of icing induced loss of control and its mitigations: AIAA-2010-8140[R]. Reston: AIAA, 2010. |
| [2] | CEBECI T.Effects of environmentally imposed roughness on airfoil performance: NASA CR 179639[R]. Washington, D.C.: NASA, 1981. |
| [3] | BRAGG M B, GREGOREK G M.Wind tunnel investigation of airfoil performance degradation due to icing: AIAA-82-0582[R]. Reston: AIAA, 1982. |
| [4] | BRAGG M B, ZAGULI R J, GREGOREK G M.Wind tunnel evaluation of airfoil performance using simulated ice shapes: NASA CR 167960[R]. Washington, D.C.: NASA, 1982. |
| [5] | BRAGG M B, COIRIER W J.Detailed measurements of the flowfield in the vicinity of an airfoil with glaze ice: AIAA-85-0409[R]. Reston: AIAA, 1985. |
| [6] | BRAGG M B, COIRIER W J.Aerodynamic measurements of an airfoil with simulated glaze ice: AIAA-86-0484[R]. Reston: AIAA, 1986. |
| [7] | BRAGG M B, SPRING S A.An experimental study of the flowfield about an airfoil with glaze ice: AIAA-87-0100[R]. Reston: AIAA, 1987. |
| [8] | BRAGG M B, KHODADOUST A.Experimental measurements in a large separation bubble due to a simulated glaze ice accretion: AIAA-88-0116[R]. Reston: AIAA, 1988. |
| [9] | BRAG M B, KHODADOUST A.Effect of simulated glaze ice on a rectangular wing: AIAA-89-0750[R]. Reston: AIAA, 1989. |
| [10] | KHODADOUST A.A flow visualization study of the leading edge separation bubble on a NACA 0012 airfoil with simulated glaze ice[D]. Columbus: Ohio State University, 1987. |
| [11] | RANAUDO R J.Performance degradation of a typical twin engine commuter type aircraft in measured natural icing conditions: AIAA-84-0179[R]. Reston: AIAA, 1984. |
| [12] | RANAUDO R J, MIKKELSEN K L, MCKNIGHT R C, et al.The measurement of aircraft performance and stability and control after flight through natural icing conditions: AIAA-86-9758[R]. Reston: AIAA, 1986. |
| [13] | POTAPCZUK M G, GERHART P M.Progress in the development of a Navier-Stokes solver for evaluation of iced airfoil performance: AIAA-85-0410[R]. Reston: AIAA, 1985. |
| [14] | POTAPCZUK M G.Navier-Stokes computations for a NACA 0012 airfoil with leading edge ice: AIAA-87-0101[R]. Reston: AIAA, 1987. |
| [15] | POTAPCZUK M G.Navier-Stokes analysis of airfoils with leading-edge ice accretions[D]. Akron: The University of Akron, 1989. |
| [16] | ZAMAN K B, POTAPCZUK M G.The low frequency oscillation in the flow over a NACA 0012 airfoil with an iced leading edge[M]//MUELLER T J.Low Reynolds number aerodynamics.Berlin: Springer-Verlag, 1989: 271-282. |
| [17] | HILTNER D W.A nonlinear aircraft simulation of ice contaminated tailplane stall[D]. Columbus: Ohio State University, 1998. |
| [18] | SIBILSKI K, LASEK M, LADYZYNSKA-KOZDRAS E, et al.Aircraft climbing flight dynamics with simulated ice accretion: AIAA-2004-4948[R]. Reston: AIAA, 2004. |
| [19] | CUNNINGHAM M A.A simplified icing model for simulation and analysis of dynamic effects[D]. Morgantown: West Virginia University, 2012. |
| [20] | 王明丰, 王立新, 黄成涛. 积冰对飞机纵向操稳特性的量化影响[J]. 北京航空航天大学学报, 2008, 34(5): 592-595. WANG M F, WANG L X, HUANG C T. Computational effects of ice accretion on aircraft longitudinal stability and control[J]. Journal of Beijing University of Aeronautics and Astronautics, 2008, 34(5): 592-595. (in Chinese) |
| [21] | 袁坤刚, 曹义华. 结冰对飞机飞行动力学特性影响的仿真研究[J]. 系统仿真学报, 2007, 19(9): 1929-1932. YUAN K G, CAO Y H. Simulation of ice effect on aircraft flight dynamics[J]. Journal of System Simulation, 2007, 19(9): 1929-1932. DOI:10.3969/j.issn.1004-731X.2007.09.009 (in Chinese) |
| [22] | 蒋天俊.结冰对飞机飞行性能影响的研究[D].南京: 南京航空航天大学, 2008. JIANG T J.Study on the effect of icing on aircraft flight performance[D]. Nanjing: Nanjing University of Aeronautics and Astronautics, 2008(in Chinese). http://cdmd.cnki.com.cn/Article/CDMD-10287-2009054739.htm |
| [23] | 王起达.结冰后飞机的纵向稳定性和操纵性研究[D].南京: 南京航空航天大学, 2009. WANG Q D.Study on longitudinal stability and maneuverability of aircraft after icing[D]. Nanjing: Nanjing University of Aeronautics and Astronautics, 2009(in Chinese). http://cdmd.cnki.com.cn/Article/CDMD-10287-2010081957.htm |
| [24] | 张智勇.结冰飞行动力学特性与包线保护控制律研究[D].南京: 南京航空航天大学, 2006. ZHANG Z Y.Research on iced aircraft flight dynamics characteristics and envelope protection control law[D]. Nanjing: Nanjing University of Aeronautics and Astronautics, 2006(in Chinese). http://cdmd.cnki.com.cn/article/cdmd-10287-2006118911.htm |
| [25] | 邵元培, 车竞, 丁娣. 大飞机机翼结冰对飞行动力学特性影响研究[J]. 飞行力学, 2018, 36(1): 12-15. SHAO Y P, CHE J, DING D. Study on the influence of wing icing on flight dynamics characteristics of large aircraft[J]. Flight Dynamics, 2018, 36(1): 12-15. (in Chinese) |
| [26] | BRAGG M B, HUTCHISON T, MERRET J.Effect of ice accretion on aircraft flight dynamics: AIAA-2000-0360[R]. Reston: AIAA, 2000. |
| [27] | SONNEVELDT L.Nonlinear F-16 model description[R]. Delft: Delft University of Technology, 2010. |
| [28] | 刘世前. 现代飞机飞行动力学与控制[M]. 上海: 上海交通大学出版社, 2014: 145-150. LIU S Q. Flight dynamics and control of modern air-crafts[M]. Shanghai: Shanghai Jiao Tong University Press, 2014: 145-150. (in Chinese) |
| [29] | US Air Force.US military specification: Flying qualities of piloted airplanes: MIL-F-8785C[S]. Philadelphia: Naval Publications and Form Center, 1982. |
