飞机机翼颤振特性测试仿真研究

于之靖 ，花 贞 ，王 烁 ，王 威 ，吴 军 ，诸葛晶昌

（中国民航大学a.电子信息与自动化学院；b.航空工程学院，天津 300300）

机翼动态变形与颤振特性是飞机性能和安全的重点考虑因素，也是气动弹性变形分析、飞行试验科目和颤振适航认证的重点测试内容之一[1-3]。因此需要获得尽可能多的机翼变形和颤振测量数据，为飞机结构设计和适航认证测试及应用提供足够的试验分析数据，以确保飞机飞行安全。

机翼颤振是机翼动态变形的极端表现形式，它是机翼在气流中受到气动力、弹性力和惯性力的耦合作用而发生的一种动态气动弹性现象，极端条件下将使飞机结构在极短的时间内遭到破坏，严重影响飞机飞行安全。从固定翼飞机问世开始，各种机翼变形和颤振问题一直备受关注[4-6]。

1 技术背景

随着对机翼变形和颤振研究的逐渐深入，相关学者提出各种关于机翼适航颤振的理论。亚临界试验技术的发展，使得各种适航测试技术得以迅速发展，完善了测试数据的来源。而颤振适航测试需要尽可能贴合真实的飞行环境来模拟测试机翼振动情况，针对性地测量激振点。颤振适航飞行试验可以准确地确定适航颤振包线，但难以控制飞行环境，而且不能综合分析颤振因素，从而激发了人们选用其他技术代替适航飞行试验。风洞试验[7]基于相对运动定理，创建地面条件模拟飞机在空中的飞行状态，但考虑到模型和环境的局限性，试验结果容易失真。梁模拟机翼试验只能用于机翼颤振特性的简单分析，精确度不高。有限元数值模拟可以计算机翼在理想或实际的环境下的各种状态，是目前常用的机翼颤振特性分析手段。偶极子格网法（DLM）计算流场相对简单，但不能模拟过于复杂的流场[8-9]。随着计算流体动力学（CFD）技术的发展，流固耦合计算方法随之得以深入研究[10-11]。

本文基于颤振适航测量的目的，采用流固耦合的计算方法分析飞机流场，进而研究机翼颤振特性，可以有效提高机翼颤振适航测量激振点选取的针对性，减少测试工作量。

2 机翼振动特性理论分析

机翼的颤振形式通常表现为扭转变形导致的结构解体、应力应变过大而导致的结构刚度不足，因此，机翼颤振的研究主要集中在其气动弹性分析以及振动响应分析。通常，机翼结构的振动测试激振区域的选取以其应力、应变或者其动态变形的集中分布区域作为依据。机翼的颤振特性可用动力学分析来确定，为完整描述飞机在飞行过程中所遇到的流固耦合[12-14]问题，研究耦合流场压力效应作用下机翼的颤振特性规律，根据Ansys的求解理论，通过机翼静态条件下的线性分析，将在耦合流场作用下求解得到的流体压力边界条件加载到机翼结构的动力学平衡方程中，得到机翼耦合条件下的动力特性方程为

其中：M、C、K 分别为质量、阻尼、刚度矩阵；Mf、R、Kf、Q分别为流体等效质量、流固耦合、流体等效刚度以及非定常气动力系数矩阵；F 结构外载荷向量分别为节点位移、节点速度、节点加速度矢量；ρf为流体密度；p流体压力矢量。

3 机翼颤振测试仿真研究方法

在上述理论分析的基础上，对B737-800机翼进行颤振特性分析。首先简化机翼模型如下：机翼材料采用LY12合金，泊松比为0.33，弹性模量为7.1×1010N/m3，密度为 2 780 kg/m3。划分网格，如图 1所示。然后分别通过Block-Lanezos算法[15]和弱耦合算法分析静态和耦合流场条件下的机翼模态，可得到对应的耦合条件下机翼的振动响应，并且通过对比可得耦合流场对机翼振动模态的影响；将耦合模态映射到其他动力特性分析中，可得到耦合流场条件下的振动特性规律，分析研究后，找到机翼颤振失稳点。

图1 B737-800有限元模型Fig.1 B737-800 finite element model

3.1 机翼耦合流场计算

借助Fluent流体分析软件，采用流固耦合的方法计算机翼流场压力。首先利用SolidWorks软件，以构建的B737-800有限元模型为基础，生成三维模拟机翼流场；采用Gambit对生成的机翼周围的流场进行体网格划分，再导入Fluent软件中。在Fluent软件中，选择基于压力的稳态求解器，SA湍流模型，选择理想气体Ideal-Gas（满足气体状态方程），并选择Sutherland定律计算其粘性；设定速度进口，自由出口为初始边界条件，选择耦合求解方法，指定对应的控制方程，获得的机翼表面的压力分布与流场计算结果如图2和图3所示。

图2 机翼压力分布图Fig.2 Wing pressure distribution

图3 流场速度分布图Fig.3 Flow field velocity distribution

3.2 机翼模态仿真研究

分别分析静态条件和耦合条件下的模态响应，观察振型变化，可得到耦合流场对模态的影响。静态条件下只考虑机翼自身的重力，采用Block-Lanezos算法分析该条件下机翼的固有振型和频率。耦合条件下，修改Fluent中获得的机翼表面压力场的数据格式，将其通过Fluent的输出接口加载到Ansys中，在Ansys中重新生成对应的机翼表面压力数据组，读取该数据组，并将其加载到相应的有限元模型表面，从而完成耦合力场的导入。首先分析加载耦合力场后的机翼静态结构力学，获得相应的预应力效应矩阵；然后进入模态求解器，采用弱耦合算法，扩展模态，分析机翼在结构刚度矩阵改变后的模态振型。图4为机翼在静态条件和耦合流场条件下的模态对比（左侧为静态条件下，右侧为耦合条件下）。

图4 静态条件和耦合条件下的前4阶机翼模态对比Fig.4 Comparison of first 4 wing modals under static condition and coupled condition

观察图4机翼在不同条件下的模态，考虑耦合力场后，机翼振型明显发生了变化：1阶模态变化微弱，但机翼的相对位移值减少了；2阶振型发生了明显的变化，在耦合流场条件下，机翼整体弯扭模态被严重削弱了；3阶模态靠近翼尖处的弯扭模态被严重削弱了；4阶振型翼尖处的弯曲模态被严重削弱了。通过对比分析得知，耦合条件下机翼的模态振幅均有所下降，但考虑到机翼的模态频率变化微弱，可以认为这种变化是由机翼在空气流场条件下的耦合作用引起的，如表1所示。

表1 两种条件下模态频率对照Tab.1 Modal frequency comparison under two conditions

本文只对这两种条件下机翼的前4阶模态进行分析研究，可发现考虑耦合力场后，机翼的模态频率略有提高，这说明在空气流场作用下，机翼产生了刚度硬化现象。两种条件下，机翼的变形位移都沿半翼展方向逐渐增大，在70%翼展位置急剧增大，且最大位移均集中在翼尖位置；弯扭振型随频率的增加变得复杂。

图5为机翼耦合流场下的各阶振型对应的应力、应变响应。可观测到，1阶振动对应的最大应力、应变集中在20%～70%翼展区，2阶振动对应的最大应力、应变集中在85%翼展区至翼尖位置；其余2阶振动对应对应的最大应力、应变集中翼尖位置。

图5 耦合条件下的应力应变云图Fig.5 Stress&strain nephogram of modal analysis

3.3 机翼谐响应仿真研究

以机翼耦合模态为基础，选取频率范围0～100 Hz，将Fluent中求解得到的耦合力转换格式，加载到机翼对应表面上，采用模态叠加法求解谐波响应特性曲线。分析频谱特性曲线，寻找振动幅度最大时的频率，作为全机翼颤振特性的重点分析频率，由此可提高颤振测试的针对性。图6为机翼X、Y、Z 3个方向的频谱图。观察其响应结果可知X、Y、Z 3个方向均是在频率为3 Hz时的机翼振动幅度达到最大，因此重点观察3 Hz时机翼的振型变化，如图7所示，分析可知在该频率下，机翼发生了严重的弯曲变形，并且最大变形位移位于翼尖处；最大应力集中于23%～70%翼展区。

图6 X、Y、Z方向频谱图Fig.6 X，Y，Z direction spectrum

图7 后处理响应云图Fig.7 Postprocessing nephogram

3.4 瞬态动力学仿真研究

图8 机翼表面压力随时间变化历程Fig.8 Change of surface pressure with time

为了获得机翼在耦合流场中的瞬态振动特性曲线以及机翼的受力分布情况，为颤振测试点的选取提供借鉴作用，通过Fluent模拟求解得到机翼表面压力载荷随时间变化的曲线，如图8所示。将该时间历程处理后导入到Ansys中，加载到相应的机翼表面，进行机翼振动瞬态动力学分析[16]，得到机翼表面压力响应，如图9所示。分析可知，机翼在0.37 s时响应最大，在该时刻对机翼模态进行扩展处理，观察最大响应时刻机翼的振动响应情况。图10为机翼在0.37s时刻的振型。可以看到，机翼在响应最大时刻发生了严重的弯扭变形，并且最大位移变形和最大应力均位于翼尖处。

图9 等效压力响应Fig.9 Equivalent pressure response

图10 0.37 s时刻机翼振型图Fig.10 Airfoil model at 0.37 s

4 实例验证

4.1 静气动弹性试验

为了预测飞机的静气动特性，保证飞行安全，中国空气动力研究与发展中心[17]对静气动弹性机翼模型进行了高速风洞试验研究。机翼模型与真实机翼几何外形、刚度均相似，采用外双梁内三梁结构，梁架与蒙皮由复合材料碳纤维加工而成。翼展长1.383 m，翼根弦长0.579 m。利用模型变形视频测量（VMD）技术对其进行变形测量，试验马赫数为0.3～1.2。机翼模型的测量点分布在20%、30%、40%、50%、60%和70%半翼展位置，如图11所示。可看出，本文的机翼颤振测试仿真结果和风洞试验的测试点分布情况基本一致。

图11 静弹性测试机翼标定点分布Fig.11 Target distribution of static aeroelastic testing model

4.2 动气动弹性试验

为了测得大展弦比机翼的动气动弹性变形和压力场变化情况，西密歇根大学的Liu等[18]提出了集成VMD系统和快速响应压力敏感涂料（PSP）以同步评估模拟机翼模型变形的方法。测试系统及机翼测试点分布示意图，如图12所示。

图12 风洞试验的视频系统及测试点分布示意图Fig.12 Video system and test point distribution in wind tunnel test

沿机翼半翼展方向分布测试点，并测量各个位置的扭转角，如图13（a）所示，曲线表明扭转角沿机翼半翼展方向逐渐变大，并且可以看出，在70%翼展处，扭转角会急剧增大。在B737-800机翼翼尖位置施加不同的扭矩，模拟仿真得到机翼扭转角沿半翼展方向的变化趋势如图13（b）所示。对比发现，仿真结果与试验结果的扭转角趋向一致。

图13 沿翼展方向的机翼扭转角分布的测量Fig.13 Spanwise wing twist distribution measurements

由机翼的静气动弹性风洞试验和动气动弹性试验可知，本文的颤振特性分析结果与机翼的测试点的分布基本一致；并且仿真的扭转角分布趋势与动气动弹性试验的扭转角分布趋势相符。

5 结语

本文针对机翼颤振适航测试的激振点区域选取的问题，以B737-800机翼有限元模型为研究对象，采用计算立体流场的方法获得机翼表面压力，将Ansys与Fluent联合起来求解机翼流场耦合-预应力耦合模态，分析机翼固有振型在耦合条件下发生的变化；再应用模态法对机翼进行谐波响应分析和瞬态动力学分析，获得机翼的频率响应和位移响应特性，进而得到机翼的振动规律，找出了机翼振动危险点，为机翼颤振适航测试提供了借鉴依据，提高了测试分析的针对性和可靠性，减少了颤振测试工作量，提高了工作效率。本文的颤振测试结果与NASA兰利研究中心进行的颤振变形测试分析中选用的激振点的分布区域一致，并且扭转角沿半翼展方向的趋向也相符，验证了机翼颤振测试仿真结果的正确性。

[1]中国民用航空局.CCAR-25-R4，运输类飞机适航标准[S].中国民用航空局,2011.

[2]赵忠良,吴军强,李 浩,等.2.4 m跨声速风洞虚拟飞行试验技术研究[J].航空学报,2016(2):504-512.

[3]孙亚军,梁 技,杨 飞,等.超临界机翼跨音速颤振风洞试验研究[J].振动与冲击,2014(4):190-194.

[4]LI X L,XU L J,TAN C,et al.Real-Time Measurement of Aerodynamic Deformation of Wing by Laser Range Finder[C]//Instrumentation and Measurement Technology Conference(I2MTC),2010:1581-158.

[5]BARROWS, D A. Videogrammetric Model Deformation Measurement Technique for Wind Tunnel Applications(Invited)[C]//AIAA-2007-1163,45th AIAA Aerospace Sciences Meeting and Exhibit, Reno, Nevada,Jan 8-11,2007.

[6]BA KUNOWICZ J,MEYER R.In-Flight Wing Deformation Measurements on a Glider[C]//RaeS Research Aircraft Opearations Conference.2014:1-12.

[7]NV NGUYEN,M TYAN,JW LEE,et al.Investigations on stability and control characteristics of a CS-VLA certified aircraft using wind tunnel test data[J].Proceeding of the Institution of Mechanical Engineers Part G Journal of Aerospace Engineering,2016,230(14):2728-2743.

[8]LIU J M,JIANG X H,QIN Y L,et al.Numerical simulation of bird impact on solid-element hollow blades by using fluid-solid coupling method[J].Journal of Aerospace Power,2010,25(10):2211-2216.

[9]宗 宁,杨广珺,吕胜利.大展弦比机翼流固耦合迭代求解模拟.航空计算技术[J].2013,43(3):13-16.

[10]朱世权,李海元,陈志华,等.弹性机翼静气动弹性数值研究[J].工程力学,2017,34(S):326-332.

[11]汪亚真,刘 明,王成军,等.可变结构小型多旋翼无人飞行器的设计[J].机械设计与制造,2016(4):186-189.

[12]仲继泽,徐自力.基于动网格降阶算法的机翼颤振边界预测[J].振动与冲击,2017(4):185-191.

[13]聂雪媛,黄程德,杨国伟.基于CFD/CSD耦合的结构几何非线性静气动弹性数值方法研究[J].振动与冲击,2016,35(8):48-53.

[14]杨秋明,朱永峰,刘 清.基于流-固耦合传热的热气防冰系统干空气飞行蒙皮温度场计算研究[J].空气动力学报,2016,34(6):721-724.

[15]DING J M,WANG Z G,ZHONG C B.Model and harmonic response analysis of composite resin concrete[J].Advanced Materials Research,2011,391-392:349-353.

[16]秦可伟,马贵春,席 园,等.基于ANSYS的机翼动力学分析[J].航空计算技术,2014(2):106-109.

[17]杨贤文,余 立,吕彬彬,等.静气动弹性模型高速风洞试验研究[J].空气动力学学报,2015,33(5):667-672.

[18]LIU T,MONTEFORT J,GREGORY J,et al.Wing Deformation Measurements from Pressure Sensitive Paint Images Using Videogrammetry[C]//AIAA Fluid Dynamics Conference and Exhibit,2011:343-349.