基于流固耦合方法的大展弦比机翼非线性颤振特性分析

马艳峰, 贺尔铭, 曾宪昂, 李俊杰, 唐长红

(1.西北工业大学 航空学院, 陕西 西安 710072; 2.中航工业第一飞机设计研究院, 陕西 西安 710089)

大型运输机、高空长航时无人机都具有展弦比大、柔性大、变形大等特点，飞机在飞行载荷作用下会产生弯曲和扭转变形，一方面结构平衡状态相对于未变形结构有一定差异，结构存在预应力，各振型和频率也有差别；另一方面，变形前后气动面的差异也会导致非定常气动力的不同。这两方面综合将使机翼的颤振特性发生变化[1]。

随着数学计算工具的进步和对复杂系统非线性空气动力学认识的加深，飞行器颤振问题得到了更深入研究。在20世纪90年代早期, Van Schoor等[2]基于完全线性理论对一个大柔性飞机进行了气动弹性分析, 讨论了结构弹性对气动弹性和飞行动力学稳定性的影响。Patil等[3]研究了高空长航时飞行器的气动稳定性和飞行力学特性, 并在考虑几何非线性的基础上, 分析了大展弦比机翼的颤振特性。Mayuresh和Dowell[4]结合理论和实验方法研究了大展弦比机翼的气动弹性相应问题, 并由此掀起了大展弦比机翼非线性气动弹性研究的热潮[5-6]。国内方面，杨智春等[7]研究了非线性大展弦比柔性机冀模型的静气动弹性特性, 并通过提取静平衡位置的剩余刚度, 计算了颤振特性。谢长川、杨超[8]基于动力学小扰动假设，建立了具有大展弦比机翼柔性飞机的全机几何非线性气动弹性稳定性分析的线性化方法和工程求解流程。崔鹏、韩景龙[9]考虑了几何非线性，进行了切尖三角翼的非线性颤振和极限环振荡的研究。

本文综合考虑机翼结构的几何非线性和气动力非线性，详细地分析了大展弦比机翼的非线性颤振特性。以平板机翼为例，计算了机翼在不同攻角下的静弹性变形、模态特性以及颤振特性。通过比较颤振计算线性解与非线性解的差别，说明对大展弦比机翼进行颤振分析时，必须同时考虑结构几何非线性和气动力非线性。

1 计算方法

常规颤振计算采用偶极子网格法计算非定常气动力，当结构发生变形，偶极子网格无法跟随结构运动形成曲面，因此它不适用于大变形的非线性颤振分析。针对大变形计算气动网格要跟随结构形成曲面的特点[10]，本文选用非定常涡格法来求解非定常气动力，涡格的形状较为灵活，沿展向2条端线也无须满足平行于来流的限制，可形成任意曲面形状[10-11]，另外一个主要差别是用涡格法计算非定常流动时采用时域法，比偶极子网格多出了从后缘到远场的尾涡网格。

1.1 涡格模型

涡格模型示意图如图1所示[12]。将拓展到远场的马蹄涡拆分成若干个独立涡格，而这些涡格一部分附着于物面上，一部分则形成尾涡，将尾涡网格布置到距翼面后缘30倍弦长为止。

图1 涡格模型示意图

设物面上的涡格共有N=Nc·Ns个(Nc为物面沿弦向的网格数,Ns为沿展向的网格数),尾涡格共M=Nw·Ns个(Nw为尾涡沿弦向的网格数),则法向洗流速度和涡强度的关系可写成如下的矩阵形式[11]:

(1)

式中：上标‘～’表示尾涡。

对于定常问题:

vn=-V∞·n

(2)

对于非定常问题:

vn=(Vs-V∞)·n

(3)

式中：Vs为物体在该控制点处的运动速度,对于非定常问题,单元法向矢量n随物体的运动而不断变化。

1.2 载荷计算

涡格上作用的空气动力与涡强度有关,它由两部分组成,一部分是与涡强度时间导数相关的Ft,另一部分是与涡强度空间导数相关的Fs。

(4)

式中：ρ∞为来流密度,S为涡格面积。Ft作用于涡格的中心,方向沿单元的法向。

当涡格位于机翼前缘时:

i=(m-1)·Nc+1,m=1,2,…,Ns

(5)

当涡格不位于机翼前缘时:

i≠(m-1)·Nc+1,m=1,2,…,Ns

(6)

式中:b为涡格沿展向的跨度,FS作用于涡格前端涡线的中点,方向沿单元的法向,通过杠杆原理将Ft和Fs分配到角点上[11]。

1.3 增量法求结构静变形

本文只考虑翼面大变形产生的几何非线性影响,不考虑材料非线性等因素,采用增量法计算,结构的增量平衡方程如下:

KTΔq=ΔP-ΔR

(7)

式中：ΔP为载荷增量,ΔR为内力增量,Δq为增量位移,KT为切线刚度矩阵:

KT=K0+Kσ

(8)

式中：K0为线弹性刚度矩阵,Kσ为几何刚度矩阵(初应力阵)。

1.4 模态分析

在变形后的平衡位置求结构响应还需求解结构在该位置的动态特性,在平衡位置的结构动力方程为:

(9)

通过模态法求解该方程,将(9)式表示为模态坐标下的减缩形式如下:

(10)

1.5 分析流程图

本文使用改进的“松耦合”方式来联合推进结构和气动方程。分析时寻求二者之间的匹配关系,颤振速压即为产生静弹变形的速压,静弹变形则为颤振分析的平衡状态。反复进行静弹和颤振计算,不断迭代使二者关系达成一致,图2给出了分析流程图。

图2 非线性颤振分析流程图

2 非定常气动力验证

建立一个展弦比为24的长直机翼,网格数为20×40,给它施加一个强迫俯仰运动:

ϑ=ϑmsin(kt)

(11)

式中:ϑm为俯仰角的振幅,k为减缩频率,t为无量纲时间,转动的轴线为机翼沿展向的中心线,减缩频率k=0.2时,升力和俯仰力矩响应历程与经典theodorson理论计算结果比较如图3所示。二者的响应曲线吻合好,验证了涡格法计算的正确性。

图3 涡格法与theodorson气动力的比较(k=0.2)

3 算例及结果分析

3.1 计算模型

计算选用的平板机翼根弦长为1 m,根梢比为1,展长为4 m,展弦比为4,前缘后掠角为8°。结构模型如图4所示,平板的厚度为0.03 mm,材料弹性模量E=71 GPa,密度为2 700 kg/m3,泊松比为0.33。主要前5阶振动频率见表1。

图4 结构有限元模型

表1 主要振型及频率

气动网格分为物面网格和尾涡网格两部分,物面网格数量为20×30,与结构网格数一致,为了便于数据交换,尾涡网格一致延伸到距后缘30倍弦长。

3.2 线性颤振分析

来流密度ρ∞=2 kg/m3,通过调整来流速度,判断结构响应是否等幅震荡来判断是否发生临界颤振,得到Vf=171.47 m/s,ωf=6.115 Hz。在颤振速度下的广义位移响应如图5所示。由图可见,颤振形态是典型的弯扭耦合型颤振,参与颤振的模态主要是机翼一弯、二弯和一扭。

图5 广义位移响应

3.3 非线性颤振分析

计算攻角取α=0°～10°,攻角增量的间隔取Δα=1°。当α1=0°时,取v1=171.47 m/s。按图2所示的流程图进行静弹和颤振的匹配计算。颤振速度和频率随攻角的变化如图6、图7和表2所示。

表2 颤振速度和频率随攻角的变化

图6 颤振速度随攻角的变化

图7 颤振频率随攻角的变化

图8 不同攻角,颤振速压下,机翼前缘的静变形

表3 不同攻角下机翼平衡位置的模态频率

如图6可见,颤振速度呈现出随攻角先略微增大,后逐渐减小的过程(拐点在α≈2°),说明结构非线性将给颤振带来不利影响。图8给出了机翼前缘线随攻角的变化,由于结构非线性的影响,随着攻角的增大,每1°攻角产生的位移增量逐渐减小。表3给出了振动模态的频率随攻角的变化,随着攻角的增大,机翼弯曲频率变化幅度不大,相比之下扭转频率减小的幅度很大,10°攻角时机翼一扭频率只有0°攻角时机翼一扭频率的61%,这也正是颤振速度随攻角减小的主要原因。

3.4 结果比较

以下讨论不考虑非线性或是部分考虑非线性对计算结果的影响。对以下3种情况的计算结果进行比较:

a) 结构非线性+气动力非线性;

b) 结构非线性+气动力线性;

c) 结构线性+气动力线性

结果如图9所示。可见,当只考虑结构非线性变形,气动网格不随平衡位置变化时,计算出来的颤振速度变化趋势与考虑气动网格变化的结果趋势完全相反,原因解释如下:

图9 考虑非线性对计算结果的影响

如图10所示,非线性模态振型基本垂直于变形后的网格,若气动网格随结构变化,则(2)式中vn=

图10 气动力线性与非线性差异示意图

|Vs·n|≈|Vs|,而当气动网格不变时,运动速度Vs和法向量n之间存在夹角,则vn=|Vs·n|<|Vs|,且随着变形的增大,vn减小的幅度越大。由(1)式可知,vn减小会直接导致涡强度的减小,相应的非定常气动力减小,从而使计算的颤振速度偏高。所以,在结构非线性求解中考虑气动面的变形是十分必要的,否则得出的结果不准确,甚至可能得出相反的结果。

4 结 论

本文针对大展弦比机翼进行了静气动弹性分析与颤振分析, 并通过考虑结构变形和内应力对系统刚度的影响, 形成了机翼几何非线性气动弹性分析的实施方法。计算结果表明：

1) 随着攻角的增大，由于结构非线性的影响，每1°攻角产生的位移增量逐渐减小。

2) 攻角的增大对机翼弯曲频率影响较小，相比之下扭转频率减小的幅度很大。

3) 颤振速度随攻角的增大呈现出先略微增大而后减小的趋势，说明结构非线性将给大展弦比机翼颤振带来不利影响。

4) 气动网格随平衡位置变化与保持不变2种情形下计算结果表明考虑气动网格变形的必要性。

从研究结果来看, 传统的线性分析方法对于大展弦比机翼不适用,因此进一步研究大展弦比机翼的几何非线性气动弹性分析方法是很有必要的。

参考文献：

[1] Dowell E. Some Recent Advances in Nonlinear Aeroelasticity: Fluid-Structure Interaction in the 21st Century[R]. AIAA-2010-3137

[2] Van Schoor M C, Von Flotow A H. Aeroelastic Characteristics of a Highly Flexible Aircraft[J]. Journal of Aircraft, 1990, 27: 901-908

[3] Patil M J. Nonlinear Aeroelastic Analysis, Flight Dynamics, and Control of a Complete Aircraft[D]. Atlanta:Georgia Institute of Technology, 1999

[4] Mayuresh J P, Dewey H H. Nonlinear Aeroelasticity and Flight Dynamics of High-Altitude Long-Endurance Aircraft[J]. Journal of Aircraft, 2001, 38: 88-94

[5] Dowell E, Edwards J, Strganac T. Nonlinear Aeroelasticity[J]. Journal of Aircraft, 2003, 40: 857-874

[6] Weihua Su, Carlos E, Cesnik S. Nonlinear Aeroelasticity of a Very Flexible Blended-Wing-Body Aircraft[J]. Journal of Aircraft,2010, 47: 1539-1553

[7] 杨智春，党会学，李毅. 大展弦比机翼非线性气动弹性特性的数值模拟研究[C]∥第十一届全国空气弹性学术交流会，昆明：2009

Yan Zhichun, Dang Huixue, Li Yi. Research of Nonlinear Aeroelastic Characteris on High-Aspect-Ratio wing[]∥CARS-2009, Kunming, 2009 (in Chinese)

[8] 谢长川，杨超. 大展弦比飞机几何非线性气动弹性稳定性的线性化方法[J]. 中国科学：技术科学, 2011, 41(3): 385-393

Xie Changchuan, Yang Chao. Linearization Method of Nonlinear Aeroelastic Stability for Complete Aircraft with High-Aspect-Ratio Wings[J]. Sci China: Tech Sci, 2011, 41(3): 385-393 (in Chinese)

[9] 崔鹏,韩景龙. 基于CFD/CSD的非线性气动弹性分析方法[J]. 航空学报, 2010, 31(3): 480-486

Cui Peng, Han Jinglong. Investigation of Nonlinear Aeroelastic Analysis Using CFD/CSD[J]. Chinese Journal of Aeronautics, 2010, 31(3): 480-486 (in Chinese)

[10] Ivan W, Chad Gibbsy S, Dowell E. Aeroelastic Analysis of a Folding Wing: Comparison of Simple and Higher Fidelity Models for a Wide Range of Fold Angles[R]. AIAA-2013-1635

[11] Bret K S, Philip S Beran. Analytical Sensitivity Analysis of an Unsteady Vortex Lattice Method for Flapping Wing Optimization[R]. AIAA-2009-2614

[12] Christopher C C, Tim F, Balakumar B. GPGPU Implementation and Benchmarking of the Unsteady Vortex Lattice Method[R]. AIAA-2013-0288, 2013