某全机跨声速颤振模型颤振特性仿真与试验验证

钱 卫，杨国伟，张桂江，郑冠男

（1.沈阳飞机设计研究所，沈阳 110035；2.中国科学院力学研究所，北京 100190）

0 引 言

对于高机动飞机，经典颤振边界在跨声速区往往存在一个凹坑，因此跨声速颤振特性是飞机设计限制的关键设计参数。根据规范［1］要求，在设计初期就应进行颤振模型风洞试验，以验证规定的速度安全余量，并验证变参数颤振的分析结果。

美国NASA兰利研究中心TDT风洞从1960年至今，对各种飞机型号进行了大量的跨声速颤振风洞试验［2］，主要解决飞机型号的跨声速颤振设计问题和跨声速颤振仿真程序的评估与验证。尤其是在F－16，B－777，F18E／F等飞机的颤振风洞试验中取得的数据，对发展颤振仿真程序作出很大贡献。科罗拉多大学对F－16的飞行状态进行了完整的CFD／CSM耦合的仿真，结构FEM模型有168799自由度，CFD网格有403919空间点，63044表面点（小攻角欧拉方程仿真）［3－6］。对于CFD跨声速颤振仿真，波音公司进一步发展了NASA基于欧拉／N－S方程的CFD程序CFL3D，进行跨声速颤振仿真，并使用HSR项目的TDT风洞模型试验结果对程序进行了修正，取得比较大的进展［7］，波音公司还使用CFL3D对B－1飞机的极限环振动（LCO）进行的全面的仿真研究［8－10］。美国在战斗机的工程颤振设计中从F－18到F－22仍然采用线化理论偶极子网格法（DLM）进行颤振计算，同时通过跨声速颤振风洞试验进行验证［11－12］。近年，ZONA公司发展了简化的欧拉方程颤振仿真软件ZEUS Code，可以进行飞机的跨声速颤振仿真［13］。

国内基于偶极子网格法（DLM）的颤振设计工作已形成比较完整的体系［14］，受跨声速风洞试验条件的限制，跨声速颤振设计的准确性受到一定的限制。近年国内在跨声速CFD颤振仿真上取得很大的进展，对全机状态，实现了跨声速颤振仿真［15］，但一直缺乏对跨声速颤振仿真结果的验证。

飞机跨声速颤振评估的难点在于强的气动力非线性造成颤振耦合机理比较复杂，目前基于非定常CFD的仿真成熟度不高，基于缩比弹性相似模型的全机跨声速颤振风洞试验技术难度很大，尤其是全机弹性结构相似模型，国内在本文之前没有先例。

针对某全机复杂耦合的颤振现象，本文成功设计了全机缩比结构相似跨声速颤振模型，在FL－26风洞开展了全机颤振风洞试验，同时使用基于N－S方程的CFD仿真程序进行了跨声速颤振仿真。通过仿真结果与跨声速风洞试验结果的对比，验证了一种全机复杂耦合的跨声速颤振形式。

1 全机跨声速颤振风洞试验

1.1 全机结构相似跨声速颤振模型

中国空气动力发展中心FL－26风洞是一座试验段尺寸为2.4m×2.4m的引射式、半回流、暂冲式跨声速增压风洞。该风洞具备一套全机颤振模型的悬浮式支持系统，是目前国内唯一具备全机颤振风洞试验条件的跨声速风洞。

根据该风洞试验段的条件，设计了全机缩比弹性相似颤振模型，模型主要使用玻璃纤维复合材料＋碳纤维复合材料＋泡沫组成主要结构，局部连接使用铝，钢等金属材料。具体见图1。

图1 全机缩比弹性相似颤振模型Fig.1 Full aircraft structuresimilar flutter model

该模型主要结构传力形式与飞机相似，利用玻璃钢、碳纤维和泡沫塑料的混合结构，实现了全机缩比弹性结构相似。蒙皮主要使用阴模，采用湿法在常温下加工分区变厚度的蒙皮构件。机翼上的梁、肋和蒙皮可以组成多闭室薄壁结构，承受机翼的载荷。梁、肋和蒙皮主要依靠粘接形成结构，见图2。

模型安装在FL－26风洞的全机悬浮式颤振风洞试验系统上，模型在风洞中的安装见图3。

图2 全机缩比弹性相似颤振模型的制作Fig.2 Manufacture of the flutter model

1.2 模型的FEM仿真与地面GVT试验

对模型进行了地面共振试验（GVT），测量模型的频率、模态，在试验中把模型安装在风洞悬浮式支持系统上，试验采用多点随机激励的方法进行，见图4。

对风洞模型进行了结构动力学限元模型（FEM）仿真，仿真模型具有约1000000自由度，见图5。模型主要频率、模态计算结果见图6。

图3 全机颤振模型在风洞中的安装Fig.3 Full aircraft flutter model in wind tunnel

图4 模型地面GVT试验Fig.4 GVT for the flutter model

图5 结构动力学仿真有限元模型（FEM）Fig.5 FEM of the flutter model

图6 主要模态计算结果Fig.6 Modes of FEM simulation

模型的主要模态试验结果与仿真计算结果振型基本一致，频率对比见表1。

表1 仿真计算与GVT结果Table 1 Result of simulation and GVT

使用亚声速偶极子网格法（DLM）进行了颤振计算，软件采用 MSC／NASTRAN－Aeroelastic，非定常气动力网格见图7。

图7 DLM法非定常气动力网格Fig.7 Mesh of DLM

2 跨声速颤振CFD仿真计算

采用CFD／CSD耦合的方法对此模型的风洞试验状态进行了仿真，非定常气动力通过求解Navier－Stokes方程（混合网格）获得，结构变形由广义结构运动方程计算。将流体方程与结构运动方程分别构造为含子迭代的二阶隐式时间推进格式，通过子迭代提高整个计算的时间精度。计算模型的气动外形和风洞模型外形完全一致。

本文采用的方法里将所有阶模态都进行参数统计，根据计算出的广义位移增长率最大的那阶模态的增长率等于0时判断为颤振临界点。而随着时间步的推进，所有模态的频率逐渐靠近为一个统一的频率，就是颤振频率。

本文采用含有附面层网格的混合网格进行计算，计算中采用的机翼表面和空间非结构网格如图8。其中模型表面的三角形网格数为555190个，附面层为9层网格，空间四面体网格数为10094920个。

为提高计算中繁重的流体／结构之间的数据转换效率，将FEM模型的结构模态数据插值到飞机物面的计算网格点上，通过结构运动方程计算物面的广义位移也转换到物面的计算网格点上。

对非结构网格，本文采用一种基于弹簧模型的动网格生成技术，通过迭代得到内网格每一网格点的变形量，将变形量叠加在原网格点的坐标值上，得到下一时间步的计算网格。

首先在马赫数为0.2时的计算颤振边界，主要与NASTRAN的DLM方法对比，主要采用变来流速度方法，计算结果见图9。然后根据风洞试验中的气流参数特性完全仿真风洞试验过程，计算马赫数取为0.9，来流速度为281.04m／s，进行变来流密度颤振计算，通过改变来流密度计算得到了振动等幅时的来流密度，见图10。

图8 模型CFD网格Fig.8 CFD mesh for the flutter model

图9 M＝0.2变来流速度的颤振计算结果Fig.9 Result of CFD simulation using variable airflow speed（M＝0.2）

图10 M＝0.9变来流密度的颤振计算结果Fig.10 Result of CFD simulation using variable airflow density（M＝0.9）

3 跨声速颤振风洞试验及综合分析

在FL－26风洞进行了全机模型颤振风洞试验，风洞试验采用固定M 数改变来流速压的吹风方式进行，通过应变测量和加速度传感器测量模型的动响应判断模型是否发生颤振，模型在风洞中的安装见图3。在马赫数0.9，速压48.1kPa，发生了颤振现象，通过风洞紧急关车使模型退出了颤振。颤振发生时的动响应见图11，所有传感器都达到同一个振动频率。

综合以上的计算和试验结果，CFD仿真在M＝0.2时与DLM法的结果颤振频率基本一致，颤振速度有一定的误差，原因可能因为DLM法对于小阻尼型颤振在阻尼为0的点误差比较大，一般取1%～3%左右的阻尼点作为颤振边界，具体对比结果见表2。

图11 M＝0.9颤振模型动响应Fig.11 Frequency response of flutter model（M＝0.9）

表2 CFD计算与DLM计算结果的对比（M＝0.2）Table 2 Comparison of CFD and DLM simulation（M＝0.2）

CFD仿真的颤振速压与风洞试验结果进行对比，在M 数0.9时误差约为9%左右，这个结果已经能够满足工程设计的需求，具体对比结果见表3。

表3 CFD计算与风洞试验结果的对比（M＝0.9）Table 3 Comparison of CFD simulation and wind tunnel test（M＝0.9）

4 结 论

通过以上的计算与试验研究获得以下结论：

（1）全机缩比弹性结构相似颤振模型，满足全机模态的模拟要求；

（2）基于Navier－Stokes方程（混合网格）的CFD仿真，在采用固定密度，M数，改变来流速度的条件下，与DLM方法对比，结果符合比较好；

（3）基于Navier－Stokes方程（混合网格）的CFD仿真，对应风洞试验状态，采用固定M数，改变来流密度的方法，与风洞试验结果符合比较好；

（4）通过全机缩比弹性相似颤振模型的跨声速颤振风洞试验以及基于Navier－Stokes方程（混合网格）的CFD仿真，验证了一种全机复杂耦合的颤振形式，在飞机非定常CFD仿真上取得了进展，可以对全机复杂耦合的颤振特性进行仿真，用于飞机的工程颤振设计。

［1］GJB 67A7－2008.军用飞机结构强度规范第7部分：气动弹性［S］.

［2］RIVERA J A，FLORANCE J R.Contribution of transonic dynamics tunnel testing to airplane flutter clearance［R］.AIAA 2000－1768，2010.

［3］ELI Livne.Future of airplane aeroelasticity［J］.Journal of Aircraft，2003，40（6）：1065－1092.

［4］FARHAT C，GEUZAINE P，BROWN G.Application of a threefield nonlinear fluid－structure formulation to the prediction of the aeroelastic parameters of an F－16fighter［J］.Computers and Fluids，2003，32：3－29.

［5］FARHAT C，PIERSON K，DEGAND C.Multidisciplinary simulation of the maneuvering of an aircraft［J］.Engineering with Computers，2001，17：16－27.

［6］GEUZAINE P，BROWN G，HARRIS C，et al.Aeroelastic dynamic analysis of a full F－16 con.guration for various flight conditions［J］.AIAA Journal，2003，41（3）：363－371.

［7］MOELJO S Hong.Evaluation of CFL3D for unsteady pressure and flutter predictions［R］.AIAA－2003－1923，2003.

［8］HARTWICH P M，DOBBS S K，ARSLAN A E，et al.Navier－Stokes computations of limit－cycle oscillations for a B－1like configuration［J］.Journal of Aircraft，2001，38（2）：239－247.

［9］ERICSSON L E.Comment on Navier－Stokes computations of limit－cycle oscillations for a B－1like configuration［J］.Journal of Aircraft，2002，39（3）：519.

［10］HARTWICH P M，DOBBS S K，ARSLAN A E，et al.Reply to comment on Navier－Stokes computations of limit－cycle oscillations for a B－1like configuration［J］.Journal of Aircraft，2002，39（3）：519，520.

［11］BRIAN Hayes W，CHARLES E.Goodman.F／A－18E／F Super Hornet flutter clearance program［R］.AIAA 2003－1940，2003.

［12］ANDERSON W D，MORTARA S.F－22 aeroelastic design and test validation［R］.AIAA 2007－1764，2007.

［13］CHEN P C，ZHANG Z C，AYAN S.Aeroelastic applications of a panel－model－based overset Euler solver to a twin－engine transport flutter model［R］.AIAA 2010－2550，2010.

［14］管德，主编.《飞机气动 弹性力学手册》［M］.航空工业出版社，1994.

［15］YANG G W.Numerical analyses of transonic flutter on an aircraft［J］.Chinese Journal of Theoretical and Applied Mechanics，2005，37（6）：769－776.