_结冰对飞机飞行动力学特性影响的仿真研究

DOI:10.16182/j.cnki.joss.2007.09.009

统 仿 真 学 报 V ol. 19 No. 9

2007年5月 Journal of System Simulation May, 2007

第19卷第9期 系

结冰对飞机飞行动力学特性影响的仿真研究

袁坤刚, 曹义华

（北京航空航天大学 航空科学与工程学院, 北京 100083）

摘 要: 采用一种计算飞机结冰后气动系数和气动导数的工程方法，引入几个结冰参数，建立结冰后飞

机飞行动力学模型。根据这些参数计算了飞机结冰后阻力系数以及典型气动导数的变化，并计算了结冰前后飞机的飞行包线；同时对结冰前后飞机的纵向、横侧动稳定性以及升降舵、副翼和方向舵阶跃操纵响应进行仿真计算。这种方法除了简单、物理概念明显、耗费少之外, 采取将外界结冰条件与飞机本体分开处理的思路, 具有广泛的适应性。结冰参数值随时间变化, 可为实时监控结冰影响提供数据。 关键词: 飞机; 结冰; 飞行性能; 飞行力学; 稳定性; 操纵性

中图分类号：V212.1 文献标识码：A 文章编号：1004-731X (2007) 09-1929-04

Simulation of Ice Effect on Aircraft Flight Dynamics

YUAN Kun-gang, CAO Yi-hua

(School of Aeronautics Science and Technology, Beijing University of Aeronautics and Astronautics, Beijing 100083, China)

Abstract: An engineering method was used to calculate aircraft aerodynamics coefficients and derivatives after icing. Several icing parameters were imported to calculate the change of drag coefficient and typical derivatives between clean aircraft and iced aircraft, and the aircraft flight envelope change of iced and un-iced situation was computed, and the small disturbance flight dynamics model was modified to study the effect of the ice accretion on the aircraft dynamics by importing the ice parameter to it, and the longitudinal and lateral stability and the elevator, aileron and rudder step response were also studied and simulated. This method is simple, and its physical concept is obvious and it needs less expense, meanwhile, the icing condition and the aircraft are considered dividually, and this idea has extensive adaptability. The icing parameter changes with time, and this can provide data for real-time calculation of ice effect.

Key words: aircraft; ice accretion; flight performance; flight dynamics; stability; control

引 言飞机结冰是飞机的表面某些部位聚集冰层的现象。当大气种含有相当数量的过冷水滴时，气流流经飞机表面，而过冷水滴惯性比空气大，不能绕过飞机，从而撞击在飞机表面，发生冻结现象，形成结冰。不同程度的结冰对飞机有不同的影响，轻则使飞机性能下降、飞行品质降级，严重时会危及飞行安全，甚至造成机毁人亡的飞行事故。结冰影响飞行性能和飞行品质的实质是改变飞机气动系数和气动导数的大小，通常是进行飞行试验得到结冰之后的气动系数和气动导数，通过MSR (Modified Stepwise Regression)方法来计算飞行性能和品质参数[1]，或者采用计算流体力学和风洞实验方法获得结冰后气动系数和导数[2]，这些方法计算量大，费用高，精度也有限。本文引入结冰参数对飞机结冰后气动系数和导数进行估算，并计算了结冰后飞行包线、动稳定性以及操纵响应的变化。

结冰条件参数(冻结系数n ，积聚参数A C ，收集效率E ) 而变化[3]。

冻结系数n 定义为: n =

m f

m imp

冻结系数的其中m f 为结冰水滴质量，m imp 为碰撞水滴质量。

值对应于该条件下结冰的类型，大的冻结系数(接近于1) 对应霜冰冰形，而小的冻结系数(接近于0) 对应光冰冰形。

积聚参数A C 为无量纲质量流量，表征在t 时刻，翼型弦长上垂直于来流方向的平面上的结冰厚度，定义为：

V (LWC ) t A C =

ρice c 其中V 为飞行速度，LWC (Liquid Water Content)为液态水含量，c 为机翼气动弦长，ρice 为结冰的密度;

收集效率E 定义为为

:

E =

m imp m p

1 方法概述

1.1 结冰后阻力系数增量计算

根据NASA 冰风洞实验结果, 阻力系数增量∆C d 随一些

其中m p 为机体投影面积上水的质量。

NASA 冰风洞实验结果表明，在冻结系数n ＝0.2的时n 偏离0.2的时候阻力系数增量下降。候，∆C d 达到最大值，

当n 从0.2增大时，结冰类型逐渐过渡为霜冰，冰形棱角逐渐消失，阻力系数增量减小，当n 从0.2开始减小时，碰撞水滴中结冰的水滴很少，阻力系数增量也减小。

阻力系数是可看作结冰时间的分段函数, 当遭遇结冰

收稿日期：2006-03-21 修回日期：2006-07-14

作者简介：袁坤刚(1982-), 男, 贵州金沙人, 博士生, 研究方向为飞行动力学与控制; 曹义华(1962-), 男, 江西都昌人, 教授, 博导, 研究方向为旋翼飞行器设计、动力学与控制等。

·1929·

2007年5月 系

统 仿 真 学 报 May, 2007

的分段线性函数g (n ) 的值，本文根据文献[3]中结冰后阻力系数的实验数据，拟合得到g (n ) 的表达式如下：

⎧22.667n ,0≤n ≤0.2

(8) g (n ) =⎨

⎩4.533−4.75(n −0.2),0.2≤n ≤1

时间t 小于10分钟时阻力系数的增量可表示为[3]

∆C d =z 1(A C E ) g (n ) (1)

式中:z 1为比例常数; g (n ) 是冻结系数n 的分段线性函数, 当n =0.2时g 取最大值，此前g 线性增加，此后g 线性减小。

当遭遇结冰时间大于10分钟时阻力系数的增量可表示为

[4]

1.5 结冰后飞行动力学模型

将结冰后的气动导数C (A ) iced 引入飞机的小扰动线化模型，就得到结冰后飞机的飞行力学模型，采用该模型进行结冰后动稳定性和动操纵性的计算。

∆C d =z 2(1−e z 3(t −10) ) +∆C d (t =10min) (2) 式中: z 2, z 3为基于在t =10min 时的∆C d 的计算系数。

1.2 结冰后最大升力系数计算

根据风洞实验结果C L max 可以表示为在某一较低迎角下升力系数/阻力系数增量的线性函数[5]:

2 仿真计算与分析

2.1 阻力系数增量计算

图1给出了在AcE =0.015和AcE =0.019两个状态下的阻力系数增量∆C d 与冻结系数n 的关系。

C L max =k 1(C L , iced −C L , clean ) =k 2(C d , iced −C d , clean ) (3)

式中: k 1, k 2为比例系数; C L , iced 为某一迎角下结冰后的升力系数; C L , clean 为未结冰飞机的升力系数。C d , iced 为某一升力系数下结冰后的阻力系数; C d , clean 为结冰前的阻力系数。

1.3 结冰后气动导数计算

根据结冰的物理机理, 飞机结冰后气动导数的变化可表示为

[3,4]

'

C (A ) iced =(1+ηice k C ) C (A ) (4) A

式中:C (A ) 为某一气动导数; C (A ) iced 为结冰后该导数值;

'

仅与气象条件有关; k C 为飞机结冰因ηice 为结冰程度参数，A

子, 它是经过修正的系数, 与飞机几何尺寸, 飞行速度，结冰条件等有关。

结冰程度参数代表了飞机遭遇到的结冰气象条件的严重程度, 采用NACA0012翼型定义为

ηice

∆C d (NACA 0012, c =3' , V =175kts , 实际结冰条件)

= ∆C d ref (NACA 0012, c =3' , 最大连续结冰条件)

图1 阻力系数增量∆C d 与冻结系数n 的关系比较

从图中可以看出，冻结系数n

线性增加，在n ＝0.2时达到最大，而n >0.2时则线性减小。和参考文献[3]结果以及实验数据相比，本文的计算结果略大，

但非常接近，这说明阻力系数增量的估算公式和本文对函数

(5)

式中分子中实际结冰条件为实际飞行中遭遇的结冰气象条件, 分母中结冰条件为FAA 附录C 里的最大连续结冰条件:液态水含量LWC =0.65g/m3; 平均水滴体积直径（Median Volumetric Diameter）MVD =20μm; 速度V =292.25km/h;结冰

'

遭遇时间t ＝10min; 大气温度T =-4℃。然而由于k C 与气象A

g (n ) 表达式的拟合是合理的。

2.2 气动导数计算

算例采用加拿大的双水獭飞机(Twin-Otter)，计算条件为：飞行高度h =1999.5m, 飞行速度V ＝70.1m/s, 结冰参数η=0.2。

在该条件下双水獭飞机结冰前后的典型气动导数如表1

表1 结冰前后双水獭飞机的无因次气动导数（η＝0.2）

气动导数 结冰前 k C A 结冰后

-5.66 -0.2809 -5.094 C z α

-19.97 -0.0676 -19.70 C zq

-0.608 -0.477 -0.55 C Z δe

-0.041 1.0976 -0.062 C X 0

-0.2943 -0.1903 -0.2831 C X α

-1.31 -0.0954 -1.18 C m α

条件和飞机构型都有关，使用ηice 在确定具体的飞机结冰后导数时存在不便。因此引入参数η和k C A ：

'

k C =A

η

k C (6) ηice

A

只与飞式中, k C A 表示飞机的某个气动导数对结冰的敏感性，机有关而与气象条件无关, 它是关于飞机翼型、展弦比、后掠角、布置形式等的函数, 对于给定飞机是定值。η的计算和ηice 相似, 只是采用实际飞机机翼气动弦长及相应的速度。它是关于大气温度、液态水含量、水滴直径、结冰遭遇时间等的函数, 只与气象条件有关与具体飞机无关。

由(4)和(6)式，可得:

C m δe

C mq

C y δr C y β

-1.74 -0.089 -1.566 -34.20 -0.1754 -33.00 0.15 -0.4 0.138 -0.6 -1.00 -0.48 -0.50 -0.50 -0.45 -0.15 -0.500 -0.135 0.015 -0.50 0.0138 0.1 -1.0 0.08 -0.18 -0.3056 -0.169 -0.12 -0.4167 -0.11

C (A ) iced =(1+ηk C A ) C (A ) (7)

C lp

C l δa

C l β C n β

1.4 函数g (n ) 的确定

在计算结冰后阻力系数的过程中，需要确定冻结系数n

·1930·

C nr C n δr

2007年5月 袁坤刚, 等：结冰对飞机飞行动力学特性影响的仿真研究 May, 2007

所示，表中还给出相应的k C A 的值。从表1可以看出，在该结冰条件下，结冰后纵向、横向的典型气动导数变化的量级并不大，几乎都在15%以内。

明显。从各个特性参数的变化来看，结冰对纵向动稳定性的影响不大。结冰后CAP 参数值增大，表明与结冰前相比，飞机操纵量不足。

结冰前后飞机横侧动稳定性特征值如表3所示。

表2 结冰前后飞机纵向动稳定性参数

特性参数 结冰前 结冰后

-0.8026±1.3302i -0.6566±1.2894i λsp

0.5166 0.4538 ξsp

1.5536 1.4470 w sp

4.7234 4.8728 T sp

0.8634 1.0555 T 1/2sp

-0.0090±0.0911i -0.0092±0.0938i λp

0.0985 0.0974 ξp

0.0916 0.0943 w p

68.9620 66.9529 T p

0.1901 0.2156 CAP

表3 结冰前后飞机横侧动稳定性特征值

特征根 螺旋模态 滚转模态 荷兰滚模态

结冰前 0.0058 -2.0319 -0.267±1.6168i 结冰后 -0.0310 -1.2240 -0.127±0.6764i

2.3 飞行包线计算

算例采用国内某型飞机，飞机质量m ＝45000kg ，机翼面积S ＝121.86m2，计算了高度从0m ～10800m 范围内，在结冰条件为LWC =0.65g/m3, n =0.2，结冰时间t ＝5min 时，飞机最小、最大平飞速度的变化，如图2所示。

1200010000H /m

[**************]0

60 100 140 180 220

V/(m/s) 图2 结冰前后的飞机飞行包线

Vmin, uniced Vmin, iced Vmax, uniced Vmax, iced

从图2可以看出，结冰后飞机的飞行包线明显向内“减缩”，随飞行高度增的增加，包线减缩的程度就越大。因此当在结冰的情况下飞行时，飞机应在结冰后的飞行包线内飞行，否则就会造成飞行事故的发生。

表3的结果表明，与结冰前相比，螺旋模态趋于稳定，而滚转模态和荷兰滚模态迅速趋于不稳定，特别是荷兰滚模态特性恶化。

2.4.2 单位阶跃操纵响应

结冰前后升降舵单位阶跃响应计算结果如图3所示。 从升降舵单位阶跃响应可以看出，结冰以后表征长周期模态的速度、俯仰角响应的稳态值略为减小，初始响应偏慢；表征短周期模态的迎角和俯仰角速度初始响应的峰值略大，响应变慢。

结冰前后副翼和方向舵单位阶跃输入响应如图4所示，图中in(1)表示副翼输入，in(2)表示方向舵输入。

副翼和方向舵阶跃输入的响应结果表明，与结冰前相比，无论是侧滑角、滚转角还是滚转角速度、偏航角速度，结冰后副翼和方向舵响应都不平稳，而且出现较大振荡，这也主要是模态特性恶化的表现。

2.4 动稳定性和操纵响应

算例采用该型飞机，计算了在飞行高度h ＝3000m ，飞行速度V =136m/s, 结冰条件LWC =0.65g/m, n =0.2时飞机的纵向、横侧动稳定性和操纵响应。 2.4.1 动稳定性计算结果

结冰前后飞机纵向动稳定性特征参数以及纵向飞行品质参数CAP （操纵期望参数）如表2所示。其中，λsp 、λp 为短、长周期模态特征值；ξsp 、ξp 为短、长周期模态阻尼比；w sp 、w p 为短、长周期模态自然频率；T sp 、T p 为短、长周期模态的周期；T 1/2sp 为短周期模态半幅时间。

从表2可以看出，与未结冰的情况相比，结冰后长周期模态的周期减小；而短周期模态的周期增大，模态特性趋于不

3

(a) 速度响应时间历程 (b) 迎角响应时间历程

·1931·

2007年5月 系

统 仿 真 学 报 May, 2007

(c) 俯仰角速度响应时间历程 (d) 俯仰角响应时间历程

图3 结冰前后升降舵单位阶跃响应

(a) 侧滑角响应时间历程 (b) 滚转角速度响应时间历程

(c) 偏航角速度时间历程 (d) 滚转角响应时间历程

图4 副翼和方向舵单位阶跃响应

3 结论

通过使用参数ηice 和η来预测不同的结冰条件对飞机性能、气动和操纵参数的影响, 并与飞行力学模型结合来研究结冰后飞机的动力学特性。计算和仿真的结果表明：

(1) 结冰参数ηice 和η可以用来计算结冰后的气动导数，并可以引入结冰后飞行力学模型进行操稳特性计算和分析；

(2) 在计算阻力系数增量的过程中，拟合而得的分段线性函数g (n ) 的表达式有较好的合理性，其精度可以满足工程需要。

(3) 结冰后飞机飞行包线向内“减缩”，飞行范围缩小；结冰后短周期模态、滚转模态和荷兰滚模态趋于不稳定，滚转模态和荷兰滚模态特性迅速恶化；纵向操纵响应不足，横侧操纵响应不稳定。

上述模型中的参数ηice 和η计算仅依据阻力增量来计算, 可以想见结冰后气动系数/导数的精度不是很高；按不同的研究目的, 可考虑增加其它相关参量来改进这一方法, 以提高动力学参量估算的精度。

参考文献:

[1] [2] [3] [4] [5]

Ratvasky T P, Ranaudo R J. Icing Effects on Aircraft Stability and

Control Determined from Flight Data [R]. NASA TM 105977, 1993. Adil Ahmed, Dafa Alla. Ice Effect on The Aerodynamic Characteristics of Airfoils [R]. AIAA paper NO. 2002-3050.

Bragg M B, Hutchison T, Merret J, et al. Effects of Ice Accretion on Aircraft Flight Dynamics [R]. AIAA paper NO.2000-0360.

Pokhariyal D, Bragg M B, Hutchison T, et al. Aircraft Flight Dynamics with Simulated Ice Accretion [R]. AIAA paper NO. 2001-0541.

Merret J M, Hossain K N, Bragg M B. Envelope Protection and Atmospheric Disturbance in Icing Encounters [R]. AIAA paper NO. 2002-0814.

·1932

·