基于非线性气动力的失速颤振计算与试验研究

戴玉婷，严 慧，王林鹏

(北京航空航天大学航空科学与工程学院，北京 100191)

自1980 年以后，科学界对动态失速和失速颤振的仿真计算有很大兴趣，各种计算模型不断发展。在过去几十年中，许多研究人员通过实验，数值模拟等方式[1 − 4]开展了研究。目前人们普遍认为没有通用的动态失速模型，即使是先进的CFD 计算[5 − 9]也不能提供一致的结果。Dimitriadis和Li[10]对NACA0012 翼型的俯仰和沉浮自由度做了一些试验研究。证明沉浮运动对失速颤振的发生没有影响，并得出失速颤振呈现出不同的分支的结论。Yabili 等[11]等利用自己制定的求解器对二元翼段低雷诺数下的动态失速和失速颤振进行了CFD 计算，并与Li 的试验结果做了对比验证。

对于工程应用而言，半经验动态失速模型在直升机和风力涡轮机空气动力学中仍然普遍存在。在这些半经验动态失速模型中，Leishman-Beddoes 模型是最受欢迎的模型之一，主要适用于0.3 Ma～0.8 Ma。Sheng 等[12]提出了一种改进的动态失速模型，该模型适用于0.3 Ma 以下气动力计算。由Galbraith 等[13]开发的一种新的失速发生标准已被用来取代基于Evans-Mort 相关性的原始标准。任勇生和刘廷瑞[14]研究了复材薄壁梁非线性失速颤振特性，得出结构阻尼在抑制颤振、增强结构稳定性方面效果明显。

失速颤振试验相关研究也从动态失速[15]、试验规律、影响因素[16]等方面有了新的进展。随着计算机，激光和相关技术的发展，二维和三维粒子图像测速(PIV)技术成功地用于测量翼型周围的流场速度[17 − 18]，分析非定常气动力产生机理。Bhat 和Govardhan[19]使用PIV 和测压元件测量了NACA 0012 俯仰振荡翼型周围流场和非定常气动力。Uruba[20]通过振荡模式分解(OPD)方法分析了俯仰振荡的NACA 0012 翼型周围的PIV速度场。Šidlof 等[21]展示了一种自振荡俯仰的翼型设计方法，该模型适用于中等雷诺数下流动引起的振动。

在Sheng 等[12]的修改模型基础上，本文对标准Leishman-Beddoes 非线性非定常气动力模型参数进行了修正，建立二维翼段气动弹性模型进行失速颤振计算，并根据计算结果设计并完成了失速颤振风洞试验，进一步验证了修正的L-B 模型可以用来进行低速翼型的失速颤振工程分析与极限环振荡评估。

1 失速颤振计算方法

1.1 非线性气动力计算方法

原始L-B 模型适用的速度范围在0.3 Ma～0.8 Ma之间，为了使该非定常气动力模型适用于低速情况，本文在Beddoes 和Sheng[3, 11]研究基础上，充分考虑翼型表面涡对气动载荷的影响，在低马赫数条件下对L-B 模型进行修正，使其适用于低速条件的非定常气动力计算。

在上升沿失速过程中，涡的形成和分离靠近前缘位置，在低马赫数时附面层在分离的时候仍有部分附着在翼型的上表面。因此，其附着效应将为上表面带来一个附加升力，称之为“超调”。超调升力是低马赫数下特有的量，该部分的升力系数可由下式计算[22]：

式中：B1=1.0为与翼型有关的参量； f为弦向分离点的位置； f′′为初始模型中对于分离点 f的一阶修正；Vx为由涡产生的法向力系数的形函数。

式中： τ为无量纲时间；Tv=7.2为沿弦向传播的时间常数；Tvl=4.5为涡传播距离的时间常数。

式中：B2=0.32为翼型相关参量； τv为涡传递的无量纲时间参数。

在下降途中，气流重新附着到翼型表面时，会受到涡的影响使得附着延迟，进而实际气动力要小于未做延迟处理的模型所得值。该差值为下降沿的超调值。αmin0=16.57°表示下降沿气流再次附着阶段的起始角度，Tr=7.02为气流再次附着阶段涡运动时间常数，当α<αmin0时气流再次进入附着阶段，翼型上表面涡运动停止，气流由分离进入完全附着阶段。在下降沿中将超调过程开始时对应的攻角设为αmin0，则超调过程与超调开始对应的攻角之间的关系见下式：

式中： Tr为重新附着过程中的时间延迟常数； q 为减缩频率。

则回程过程中，在下降沿中法向力系数的超调值可写成：

式中， Vxr为 回程中基于翼型的法向力系数的形函数，其计算式为：

式中， τr为回程的无量纲时间。

对应的力矩系数：

低速修正模型最终的法向力系数和力矩系数：

1.2 失速颤振计算方法

图1 为二维翼段气动弹性系统，该系统具有俯仰和沉浮自由度。综合1.1 节非定常气动力计算方法，则该系统的气动弹性方程可以写为：

式中： θ为翼型的俯仰角； h为翼型的垂向位移；m 为翼型单位长度质量； Kh和 Kθ分别为沉浮和俯仰弹簧刚度； Ch和 Cθ分别为沉浮和俯仰阻尼系数；Sθ为质量静距； Iθ为翼型相对1/4 弦线处的极惯性矩； Qh和 Qθ分别为外部气动力和外部气动力矩；其表达式在1.1 节中已求得，具体如下：

图1 二元翼段气动弹性系统Fig.1 Two-dimensional airfoil aeroelastic system

式中：CL和CM分别为升力系数和气动力矩系数；c 为弦长；l 为展长；xh为弹性轴距1/4 弦线的距离；ρ 为空气密度；V 为来流速度。

L-B 模型中的攻角α 和俯仰减缩角速率q 可由下式进行计算：

将式(3)改写为时域状态空间方程形式，可得到：

式中：

当给定初始条件x ，利用显式或者隐式积分可由求得其时域响应：在本文中，使用四阶定步长龙格库塔法求解该状态空间方程。

1.3 数值计算及结果分析

为验证低马赫数下的修正模型，仍然选用NACA0012 翼型。并将计算结果与原始模型和试验以及Sheng 的模型的结果[3]做了对比。验证中来流速度为0.1 Ma。对比图中给出了深度失速条件(攻角α=15◦+10◦sinωt)下修正L-B 模型的计算结果与试验对比，对比结果如图2 所示。

图2 攻角为α=15◦+10◦sinωt 时翼型的法向力与力矩系数Fig.2 Normal force and moment coefficient of airfoil when α=15◦+10◦sinωt

从图2 可以看出，在低马赫数下，原始L-B模型计算数据较试验数据有较大差异，主要体现在超调位置，原始模型和试验数据重合性较差。而修正之后的L-B 模型在超调位置上与试验数据吻合的较好。说明修正之后的L-B 模型有较好的精度，能够用来计算翼段在低速情况下的气动力特性。模型中俯仰力矩系数还受 B2、 Tvl等因素影响，所以力矩对的不是特别好。

1.4 失速颤振结果分析

本文选用参考文献[10]中的试验条件和试验结果作为数值算例的参数，来验证失速颤振计算结果的正确性。文献展示了在低马赫数下基于NACA0012 翼型的失速振荡试验。该试验的重要参数见下表1。

表1 试验重要参数Table1 Important parameters of the test

本文计算了二维翼段失速颤振的俯仰分岔图和极限环振荡曲线，并与试验[10]的计算结果进行了对比，对比如图3 所示。

图3 气动弹性系统的俯仰分岔图和极限环振荡曲线图Fig.3 Pitch bifurcation diagram and limit cycle oscillation graph of aeroelastic system

从分岔图中可以看出，当风速为12.8 m/s，翼段开始周期性俯仰振荡，即开始分岔。当风速较小时，尤其当风速V 低于15 m/s 时，计算数据与试验数据吻合较好，但是在风速V 大于15 m/s 时，俯仰角振荡幅值明显高于试验值。且从试验数据来看，随着风速的增加俯仰振荡幅值未有明显变化。这与结构中沉浮刚度和俯仰刚度可能存在非线性或阻尼非线性等因素有关。

2 失速颤振风洞实验

2.1 试验设计

为了探究失速颤振现象的形成机理，根据风洞的参数，设计并制作了二维翼段失速颤振的风洞试验模型，进行了多种状态下的风洞试验。失速颤振试验在北京航空航天大学D4 开口风洞中完成，开口段尺寸为1.5 m。开口风速范围小于70 m/s。

翼段翼型采用与仿真条件相同的NACA0012翼型，弦长为0.36 m，翼展为1 m。翼段翼肋之间相距0.125 m。另外，在翼弦向20%和70%处分别布置前梁和后梁。厚度和展长能保证最大瞬时风洞堵塞因子小于5%，满足试验要求。在前梁两端处安装加速度传感器，转轴处安装角位移计。

翼段质量为2.363 kg，转轴位于弦长的37.4%位置处，重心处于弦长的39.04%位置，绕转轴的转动惯量 Iθ为0.135 kg·m2。支撑系统为转轴滑轨系统。图4 为试验装置模型，其中，图4(a)为试验系统的原理图，图4(b)为试验系统的实物图。图4(a)中：1 为八个线性弹簧；2 为四只初始攻角调节器；3 为两组滑块系统；4 为沉浮刚度调节的线性弹簧；5 为支撑系统；6 为两组转轴摇杆机构；7 为角位移传感器；8 为翼段转轴；9 为加速度传感器；10 为翼段。最终翼段沉浮刚度 kh为4896.595 kg/m，俯仰刚度 kθ为135.28 N·m/rad。角位移传感器的精度为0.1 度。

试验目的是通过给定初始攻角，测量翼段发生颤振的速度，以及发生颤振后翼段的俯仰振幅。在试验过程中，通过调节拉杆两端连接器(图4(a)中2)的长度，由角度传感器测量结果，将初始攻角调节到要求的位置。然后给定一系列稳定风速，测量该初始攻角和给定速度下的角加速度和角位移，判断失速颤振特性。

2.2 试验结果与仿真对比分析

试验就四个初始攻角进行了测试，依次为7.5 度、10 度、13 度、16 度。试验风速从6 m/s～20 m/s。图5 为16 度初始攻角时，在翼段转轴处的角度传感器在不同风速下测得角位移动态曲线图。从下列图中可以看出在风速为6 m/s 时，曲线几乎平直，幅值近乎为0。风速增加时，角度幅值信号随风速的增加而增加，且呈等幅震荡形式。将不同风速下的俯仰角振荡曲线幅值提取出来，则可生成攻角为16 度时翼段俯仰角度幅值随风速变化的分岔图，分岔图见图6。图6 为不同初始攻角下翼段俯仰角随风速变化的分岔曲线图，从图中可以看出，在初始攻角为16 度时，翼段发生俯仰自由度等幅振荡的风速为6.5 m/s，即分岔风速为6.5 m/s。在初始攻角一定的情况下，俯仰角的幅值随风速的增加而增加。另外，当初始攻角增大时，俯仰角的分岔速度明显降低。

图4 翼段气弹系统风洞试验模型Fig.4 Wind tunnel test model for wing airframe system

根据第2 节二维气动弹性计算模型与试验模型参数，对试验模型失速颤振特性进行了数值模拟。模型的阻尼大小按照经验取0.1 倍的刚度。图7 为初始攻角为16 度时风洞试验与计算模型在不同风速下的对比图。其中：图7(a)为风速8 m/s时俯仰角响应；图7(b)为风速10 m/s 时的俯仰角响应。从图中可以看出，仿真计算振幅和试风洞试验振幅吻合度较好。在风速较小时，受风速调节控制系统、开口外界环境影响和支撑支座的影响，气动弹性时域响应有一定的波动。图8 为初始攻角16 度、风速10 m/s 试验和仿真的振荡频率曲线对比，从图中可以看出，振荡频率吻合较好。图9 为不同初始攻角下，试验测得的翼段俯仰角等幅振荡的幅值随风速变化的分岔曲线与数值计算的分岔曲线的对比图。从图中可以看出，除了在初始攻角为7.5 度时稍有偏差，数值计算的结果和试验测得结果吻合较好。尤其是初始攻角较大时，吻合度较高。说明修正后的L-B 模型能够较好应用于大攻角失速颤振分析。

图5 初始攻角16 度不同风速下翼段失速颤振试验曲线Fig.5 Wing stall flutter test curve of different wind speed at initial angles of 16 degrees

图6 不同初始攻角下翼段俯仰角随风速变化的分岔曲线Fig.6 Bifurcation curve of wing pitch angle with wind speed under different initial angles of attack

图7 初始攻角16 度试验和数值模拟计算曲线对比Fig.7 Test and numerical simulation calculation comparison curve at initial attack angle of 16 degrees

图8 初始攻角16 度风速10 m/s 试验和数值振荡频率曲线对比Fig.8 Comparison curve of test and numerical oscillation frequency with initial attack angle of 16 degrees and wind speed of 10 m/s

图9 不同初始攻角下翼段俯仰角随风速变化的试验和计算分岔曲线对比Fig.9 Comparison of test and calculated bifurcation curves of wing pitch angle with wind speed under different initial angles of attack

3 结论

本文着眼于非定常气动力的不确定性对气动弹性带来的影响，对标准的气动力动态失速模型进行了低速修正，建立了用于考察失速颤振现象的二维气动弹性模型，根据风洞试验的结果，可得出以下结论:

(1) 修正的L-B 模型可以用于二元翼段或低速大展弦比飞机大迎角失速时的非定常气动力分析。修正模型在低速状态下和试验数据吻合较好，尤其是在去程和回程的超调部位，明显优于原始模型。

(2) 设计了二元翼段大攻角失速颤振试验，试验表明初始攻角越大，失速颤振分岔速度越低，极限环振荡幅值越大。风洞试验结果与非线性失速颤振理论计算结果吻合较好。说明基于修正L-B失速气动力模型具有工程上可以接受的精度。