非对称结冰情形下飞机横航向系统稳定性分析方法研究

郭林亮, 王国智, 裴彬彬, 祝明红, 代胜吉

(1.中国空气动力研究发展中心低速空气动力研究所，四川绵阳，621000; 2.空军工程大学航空工程学院，西安，710038；3.32035部队，西安，710000)

飞机结冰会破坏飞机外部流场，导致飞机的空气动力学特性恶化，破坏飞机的稳定性和操纵性，驾驶员如果缺乏对结冰风险态势的准确感知，操纵失效会诱发飞行事故。为了实现飞机结冰情形下的安全飞行，各国学者针对结冰后飞机气动特性、结冰致灾物理链路、容冰安全飞行以及结冰情形下的风险评估等问题展开了大量的研究。NASA借助于风洞试验和空气动力学仿真计算等手段，分析研究了飞机结冰部位(机翼、机身、平尾等)、冰型(明冰、霜冰、混合冰)以及结冰程度(轻度、中度、重度)等因素对飞机动力学特性的影响[1-3]。文献[4]提出飞机结冰多重安全边界的概念，梳理了飞机结冰致灾物理链路的关键因素，分析了结冰对飞机飞行品质及操纵特性的影响。文献[5～8]基于非线性的结冰情形下复杂系统模型，应用微分动力学稳定域理论，分析结冰对飞机的动力学安全边界的影响。文献[9]构建了结冰情形下人-机-环实时仿真系统，致力于为结冰情形下的飞行事故重现和飞行控制律验证提供平台支撑。文献[10～12]基于控制律重构方法，设计了飞行安全边界保护系统，实现容冰飞行。文献[13]提出安全窗概念，通过复杂动力学仿真的方法获得结冰情形下安全操纵指引，提升结冰情形下的飞行员情境感知意识，辅助飞行员制定安全操纵策略。文献[14]结合可达集和极值理论，评估飞机着陆阶段的风险概率。文献[15～16]基于Copula理论，描述风险极值参数的相关性，评估带冰飞行的事故风险概率。

上述研究主要是针对对称结冰情形下飞机纵向通道的安全保障问题。但是飞行过程中如果出现除冰系统部分失效导致的非对称结冰问题，不仅会使得飞机升力减小，阻力增大，而且会诱导产生附加的滚转力矩和偏航力矩，影响飞机横航向系统稳定性，严重时导致飞机进入失控的危险状态。因此，本文搭建了非对称结冰情形下的飞机横航向系统模型，运用蒙特卡洛算法划分计算状态点集，通过动力学仿真计算的方式确定系统的稳定状态点集。根据Lyapunov稳定性理论，系统中各稳定状态点的敛散性是相互独立的，通过并行计算的方式，运用连续推进算法，计算稳定状态点附近的稳定域边界，从而构建多参数耦合的系统稳定域。该稳定域可为非对称结冰情形下驾驶员风险态势感知能力提升以及安全操纵提供理论参考。

1 非对称结冰情形下飞机横航向系统模型

非对称结冰情形下的飞机横航向系统模型如图1所示。模型主要由4部分组成：

图1 非对称结冰情形下飞机横航向系统模型

2)舵机模型。模拟舵机在控制指令作用下，输出的实际偏转量；

3)横航向增稳控制律。保证飞机横航向姿态稳定与控制，实现期望的动态响应；

4)非对称结冰影响模型。模拟结冰对飞机气动力的影响。

1.1 飞机横航向运动学模型

根据飞机横航向运动学方程组[17]，飞机横航向运动状态方程如式(1)所示：

(1)

式中：α为迎角；β为侧滑角；p为滚转角速度；r为偏航角速度；ψ为滚转角；δα为副翼偏角；δr为方向舵偏角；V为空速；Yi、Li、Ni、Yδi、Lδi、Nδi均为横航向运动大导数，如表1所示。

表1 横侧向运动大导数

1.2 舵机模型

图2为舵机模型。

图2 舵机模型

横航向姿态控制主要是通过作动器控制副翼和方向舵实现，其动力学模型如式(2)所示：

(2)

式中：T为作动器时间常数。副翼的作动器响应时间常数设定为Ta=40 s，方向舵的作动器响应时间常数设定为Tr=35 s。

副翼和方向舵的作动器工作幅度和工作速率受到实际因素的限制，根据经验给定副翼和方向舵作动器的工作幅值及速率限制参数，如表2所示。

表2 舵机参数

1.3 非对称结冰影响模型

通过结冰影响因子[18]表征冰型对飞机气动力导数以及操纵导数的影响，如式(3)所示：

CA(ice)=CA+ηki,iceCA

(3)

式中:CA(ice)为结冰后气动导数；CA为干净外形的气动导数；η∈[0,3]为结冰影响因子，表征结冰严重程度，η=0表征无结冰，η=1表征轻度结冰，η=2表征中度结冰，η=3表征重度结冰；ki,ice为敏感系数，如表3所示，表征气动参数对结冰的敏感特性，其值与结冰程度以及冰型无关，与机翼相对厚度、展弦比、后掠角等相关，可通过结冰风洞试验或试飞试验获得。

表3 飞机气动参数的结冰敏感系数

本文研究的非对称结冰情形为：右侧机翼除冰系统故障，机翼左右两侧冰型出现非对称特征。受机翼左右两侧冰型的影响，飞机机翼两侧的流场出现较大差异，两侧机翼的升力及阻力会出现较大的差值，从而产生附加的滚转力矩和偏航力矩，如式(4)所示。

(4)

式中：dmgc为机翼平均气动力弦线到机身中心线的距离。

根据式(4)计算，附加的滚转力矩系数ΔClice和偏航力矩系数ΔCnice如式(5)所示：

(5)

则非对称结冰情形下的滚转合力矩系数和航向合力矩系数如式(6)所示：

(6)

将附加滚转力矩和偏航力矩以及结冰情形下的气动参数带入飞机横航向运动学模型，即可进行非对称结冰情形下的飞行仿真。

1.4 横航向增稳控制律

图3 横航向增稳控制系统

2 系统稳定性分析方法

根据系统理论，稳定性表征研究对象受扰偏离后自主恢复初始平衡状态的能力。系统稳定性分析方法以Lyapunov稳定性理论为基础，相继出现了平方和法、分叉分析法[19-20]、相平面法[21]、流形法[22]、正规形法[23]以及可达集法[24]等，这些研究方法更加关注局部的稳定性分析，缺乏对系统参数耦合条件下的系统整体稳定性分析及描述。本文运用蒙特卡洛算法划分计算状态点集，并通过动力学仿真计算的方式确定系统的稳定状态点集，从状态点敛散性判断的角度出发，运用连续推进算法，计算稳定状态点附近的稳定域边界，从而构建多参数耦合的系统稳定域。

2.1 稳定性基本理论

自治系统如式(7)所示：

(7)

式中:x为n维状态向量；f(x)为n维的状态方程。

‖x-xe‖≤ε

(8)

2.2 基于连续推进算法的稳定域构建

基于连续推进算法的稳定域构建流程如图4所示。

图4 连续推进算法流程

算法流程如下：

1)用蒙特卡洛算法划分计算状态点集X={x1,x2,…,xm},m∈R，通过动力学仿真手段获取稳定状态点集Xe={xe,1,xe,2,…,xe,m}；

5)重复步骤3)和步骤4)至∀li

精度lacc=0.01条件下，算法效能对比如表4所示，与传统的计算方法相比较，本文提出的连续推进算法，在稳定域求解方面具有精度可调、计算效率高的突出优点，原因如下：

表4 算法效能对比

1)计算终止条件收敛于稳定平衡点变化为收敛与稳定邻域，避免了不必要的耗时；

3)采用并行计算方法实现计算效能的提升。

3 案例仿真及分析

本文用于案例仿真分析的背景飞机基本参数如表5所示。控制律参数取值如表6所示。

表5 背景飞机的基本参数

表6 控制律参数

仿真案例场景设定为：飞行高度H=3 000 m，飞行马赫数Ma=0.4，右侧机翼除冰系统故障导致飞机右侧机翼重度结冰，结冰因子η=3。

稳定域计算范围如式(9)所示：

(9)

基于蒙特卡洛算法，以偏航角速度r∈[-0.4∶0.01∶0.4]、滚转角速度p∈[-0.6∶0.01∶0.6]、以及侧滑角β∈[-0.6∶0.01∶0.6]将稳定域计算范围划分成80×120×120个计算状态点。通过复杂动力学仿真获取无冰情形下的稳定状态点集见图5，非对称重度结冰情形下的稳定状态点集见图6。

图5 无冰情形下的的稳定状态点集

图6 非对称重度结冰情形下的稳定状态点集

基于连续推进算法获得了机翼非对称重度结冰与无冰情形的系统稳定域包络，见图7。对比分析机翼非对称重度结冰前后的系统稳定域可知：机翼非对称结冰情形会导致飞机横航向系统稳定域偏离并萎缩，机翼非对称重度结冰情形下的稳定域萎缩78.65%，稳定域中心点由(0,0,0)偏移至(0.005 1,0.054 5,-0.003 1)。

图7 机翼非对称结冰与无冰情形的稳定包络对比

截取r=-0.1 rad/s，r=0 rad/s，r=0.1 rad/s的稳定域，获得填充的稳定域包线，如图8～10所示，红色填充面表征机翼非对称结冰情形下的飞机横航向系统稳定域，绿色填充面表征干净构型的飞机横航向系统稳定域。据图分析对比可知：机翼非对称结冰不仅导致系统稳定域大幅缩减，同时导致横航向系统稳定域参数耦合关系加强。填充的稳定域面积对比见表7。

图8 r=-0.1 rad/s稳定包线

图9 r=0 rad/s稳定包线

图10 r=0.1 rad/s稳定包线

表7 稳定域面积

4 结语

本文搭建了飞机非对称结冰情形下的横航向系统模型，基于Lyapunov稳定性理论，采用连续推进算法，获得多参数耦合的横航向系统稳定域，可以直观地反映机翼非对称结冰对飞机横航向系统稳定性的影响。仿真分析结果表明：机翼非对称结冰导致了横航向系统稳定域明显缩减，稳定域参数的耦合特征增强，警示非对称结冰情形下飞行员干预操纵时，避免粗暴操纵飞机舵面，防止飞机状态跃出稳定域边界进入失控的事故状态。