基于CFD/CSD方法的跨声速静气动弹性数值模拟应用研究

郭洪涛, 陈德华, 张昌荣, 吕彬彬, 王晓冰(. 中国空气动力研究与发展中心 空气动力学国家重点实验室, 四川 绵阳 6000；. 中国空气动力研究与发展中心 高速空气动力研究所, 四川 北川 6763)

0 引 言

现代先进大型飞机的巡航速域通常在跨声速范围。飞机升力线斜率在跨声速时最大且变化最剧烈，此时流固耦合作用也最敏感，结构上即使一个细微变化也可能对飞机的气动特性造成很大影响[1]。试验表明：在相同速压情况下，与亚、超声速比较，跨声速范围内静气动弹性造成的飞机升力损失最大，使焦点前移最甚[2]。随着气动弹性计算技术的进步，近年来国内外在这方面已经开展了大量的研究工作。

70年代，基于线化位流理论发展起来的面元法开始广泛应用于飞机的气动弹性计算，现在采用该方法的NASTRAN等商业软件仍被许多工程单位广泛使用[3]。1974年，Ballhaus与Steger提出了一种二维全隐式有限差分格式构造方法，这个方法导致了第一个专业气动弹性计算软件LTRANS的诞生。80年代初期，跨声速小扰动方程开始被用于跨声速气动弹性计算问题，波音公司开发了XTRAN3S，实际上是把LTRANS从二维扩展到了三维。整个80年代，全速势方程加边界层修正方法成为民用飞机设计的主流CFD方法。1985年，为了满足未来空军战斗机的需求，AFWAL与NASA合作开发了ATRAN3S，主要是考虑到小展弦比战斗机机翼的气动弹性计算[4]。上述方法对气动弹性计算的发展做出了许多开创性的贡献，但不足的是它们都是基于位流理论假设。当遭遇比较复杂的流动现象，比如当有较强的激波通过翼面或者有流动分离时，这种假设就不太适合了。80年代到90年代，基于Euler方程的各种激波捕捉格式得到发展，在被应用到飞行器气动力计算的同时，也被运用到飞行器气动弹性计算中来[5-6]。随着90年代基于雷诺平均N-S方程(简称RANS方程)求解方法的成熟，既可以考虑跨声速流动中的激波间断又可以考虑气动力湍流黏性效应的气动弹性数值计算方法也得到了有效的发展[7-8]。Robinson与Bantina在CFL3D中加入了气动弹性计算的相关模块[9]，Ames研究中心也开发了专门的气动弹性计算软件ENSAERO，这些软件采用Euler/N-S方程求解气动力，采用隐式B-W格式，B-L、k-ε等湍流模型，采用模态叠加法求解结构运动方程，可用于解决复杂外形的气动力计算和气动弹性分析研究[10]。但是，这类方法基于模态运动方程求解结构变形，对于静气动弹性计算而言效率太低。近年来，具有更高精度的大涡模拟方法和直接数值模拟方法在国内外也被广泛研究[11-15]，但这两种方法因为计算量巨大，目前仅适合于一些基础理论研究，并不适合于飞机等具有复杂边界条件物体的静气动弹性数值计算。

综上所述，对于跨声速流动且物面边界条件复杂的情况，目前还没有一套权威、统一的计算分析方法和代码，还值得深入细致地研究。本文通过耦合求解RANS方程与静气动弹性平衡方程，采用结构化动网格技术和多物理场数据插值技术，实现了跨声速复杂外形飞行器的静气动弹性数值预测分析，并基于典型风洞试验结果验证了计算方法、代码的有效性。此外，基于数值模拟结果，分析了静气动弹性对典型大型飞机机翼几何变形特性、表面压力分布以及气动性能的影响。

1 计算方法

1.1 CFD、CSD控制方程

控制方程采用时间相关的三维守恒型可压缩RANS方程，在一般曲线坐标系(ξ、η、ζ)下，其无量纲形式为：

(1)

式中：t为时间，Q为守恒变量，F、G、H为无黏矢通量，Fv、Gv、Hv为黏性矢通量。

湍流模拟选用一方程的S-A模型。物面为绝热无滑移边界条件，远场为压力远场无反射边界条件。为了加速CFD计算收敛，应用了多重网格技术。

物面的结构弹性变形采用静气动弹性平衡方程来求解。该方程为常系数齐次线性方程组，实际上相当于模态运动微分方程组没有了时间导数项，相对而言，计算量大大减小。即：

us=CFs

(2)

式中：us为结构点的变形位移矢量，C为结构点柔度影响系数矩阵，Fs为作用在结构点上的气动力、质量力、发动机推力等合力矢量。

1.2 流固耦合数据插值方法

大型飞机静气动弹性属于三维空间的小变形情况，综合内存占用率、计算效率以及插值精度考虑，采用三维薄板样条插值方法(TPS)是比较适宜的，TPS的详细计算方法可见参考文献[16]。

1.3 动网格生成方法

针对结构化流场计算网格，本文采用RBF结合TFI的方法来进行动网格生成。RBF可以视为面样条函数插值方法的三维扩充。其插值公式为：

(3)

其中，ri=(xi,yi,zi)为位移已知点，其数目为n，φ为关于空间距离‖r-ri‖的基函数，本文取φ(‖r-ri‖)=‖r-ri‖3、ψ=b0+b1x+b2y+b3z，插值公式的系数可通过已知点ri的位移di和平衡条件得到。详细计算方法见参考文献[17]。

基于RBF方法插值得到各网格块边内点的位移后，各网格块面内点的位移可基于网格块边通过TFI方法插值得到。同理，各网格块内部点的位移也可基于网格块面通过TFI方法插值得到[18]。

2 计算方法验证

2.1 计算模型与CFD网格

本文对某大型飞机翼身组合体模型进行了静气动弹性数值计算，并与风洞试验结果进行了对比验证。图1给出了计算模型的外形轮廓及CFD物面网格。该模型采用了超临界扩张后缘翼型以及大展弦比后掠下单翼布局，忽略了短舱/挂架、襟翼滑轨整流罩等部件，计算网格单元规模为500万。

(a) 外形轮廓

(b) 初始物面网格

2.2 计算方法验证

验证计算状态为：Ma=0.78，q=35 kPa，雷诺数Re=7.3×106，迎角α为-4°～8°。

静气动弹性计算是一个CFD/CSD的反复耦合迭代计算过程。对于工程应用来说，需要迭代计算效率高，收敛速度快。基于本文方法，一个计算状态只需要6～7个变形迭代步就基本上收敛了，如图2所示。在计算过程中，由于流场计算网格的拓扑结构与网格数量不会改变，于是每个变形步后的CFD计算都可以利用上个变形步的CFD结果作为初值进行续算，因此，中间变形迭代步的CFD计算只需要较少的计算步数。一般情况下，本文CFD/CSD方法的计算量是对应的刚性模型CFD计算量的1.5倍～2倍。

图3给出了本文计算结果与相应风洞试验结果的机翼升力系数CL和弯曲变形Δz的对比结果。

图2 静气动弹性变形迭代计算收敛历程Fig.2 Iteration convergence history of static aeroelastic deformation

(a) CL～α曲线

(b) 不同迎角下Δz～y/b曲线

图3中给出的CL进行了计算与试验的相关性修正。主要是修正了边界层转捩模拟不一致、风洞试验洞壁干扰、计算模型与试验模型外形差异等，但受限于当前的技术手段，还无法对CFD中的湍流模拟不足与分离特性模拟不准以及实际机翼变形可能存在的非线性因素等影响进行修正。总体上看，本文计算结果与风洞试验结果的一致性较好。CL的差异主要出现在翼面气流出现分离之后的CFD模拟误差所致；Δz的差异主要存在于小变形范围，根据机翼变形应该是连续一致的规律性推测，应该是小变形时测量设备的测量误差相对较大。

3 跨声速静气动弹性数值模拟分析

大型飞机为了提高气动效率，其巡航Ma通常位于跨声速范围。在跨声速时，即使机翼迎角很小，由于激波与边界层的强相互作用，黏性影响也不应该被忽略。因此，必须利用N-S方程计算气动力。不失一般性，本文重点考察了Ma=0.78时算例模型的静气动弹性特性。

3.1 几何变形影响

图4给出了Ma=0.78、q=35 kPa、α∈[-4°～8°]时机翼的Δz与Δε(弹性扭转角)沿展向变化曲线。

(a) 机翼弯曲变形

(b) 机翼弹性扭转变形

可以看出，当迎角为正时，静气动弹性使机翼产生向上的弯曲挠度变形，并在机翼的顺气流剖面产生负的弹性扭转角。从图中还可以看出，内翼段的变形量较小，在机翼拐折处以外机翼的弹性变形量越来越大，到翼尖处变形量达到最大值。这种变化规律与机翼刚度沿展向逐步减小的分布特性是相符的。根据参考文献[19]可知，对于后掠机翼来说，这种变形特征将减小机翼沿展向各剖面的当地迎角，从而影响载荷分布，改变机翼气动特性，产生所谓的静气动弹性效应。例如，依据本文算例，在Ma=0.78、q=35 kPa、α=2°时，翼尖的弯曲变形可以达到机翼半展长的4%，翼尖剖面的顺流向弹性扭转角与顺流向迎角变化量分别可达-2°与-4°。

3.2 压力分布影响

图5与图6分别给出了Ma=0.78、q=35 kPa、α=2°时刚性与弹性机翼的表面压力等值线分布图及不同展向位置的压力系数分布图。从图6可以看出，与机翼的几何变形特性相对应，静气动弹性对翼根的压力分布特性影响较小，对翼尖的压力分布特性影响较大，比较明显的影响是减少了靠近外侧机翼上表面的负压范围。

图5 静气动弹性对机翼表面压力的影响Fig.5 Effects of static aeroelasticity on wing surface pressure

(a) y/b=20%

(b) y/b=95%

从图6中机翼不同剖面变形前后的压力系数对比可知，越靠近翼尖时压力系数分布变化越大。这是由于在机翼小迎角(表面未出现大面积的气流分离)情况下，后掠弹性机翼的气动载荷使得机翼向上弯曲，同时由于负的弹性扭转角，导致顺气流方向的当地迎角减小，使得弹性变形后机翼的前缘吸力峰值降低，同时使得弹性机翼的升力系数小于刚性机翼的升力系数。依据本文算例，在Ma=0.78、q=35 kPa、α=2°时，弹性相对于刚性机翼的CL降低了15%以上。因此，在跨声速时静气动弹性对大展弦比后掠翼的气动载荷影响不应该被忽略。

3.3 气动性能影响

图7给出了Ma=0.78时弹性机翼的典型静气动弹性影响因子随速压变化曲线。可以看出，速压对弹性机翼升阻特性及纵向静稳定性都有不同程度的影响。其中，K1(升力线斜率影响因子)随速压基本上呈线性减小趋势，在q=35 kPa(对应巡航飞行高度)时影响因子大约为0.94，但这种线性规律性将在速压达到某一较大值后发生变化，当q>65 kPa时，K1～q曲线开始出现拐折。K2(纵向静稳定性裕度影响因子)随速压增大总体上呈现出非线性减小趋势，也即是说，速压越大，弹性机翼的静稳定性裕度越小。但是，在q>65 kPa后，K2～q曲线开始出现上翘，即影响特性发生了逆向变化。对于K3(巡航升阻比影响因子)来说，尽管静气动弹性会降低机翼的升力，但是由于降低升力的同时也降低了升致阻力，因此实际上会增大升阻比。从图中可以看出，K3在小速压情况下呈线性变化规律，但在q>45 kPa以后非线性特征逐渐显现。

图7 静气动弹性对机翼气动性能的影响Fig.7 Effects of static aeroelasticity on wing’s aerodynamic performances

综合跨声速静气动弹性对大型飞机纵向气动特性的影响来看，当速压较小时，静气动弹性影响量随速压呈线性变化，当速压增大到一定程度后，变化规律将出现非线性特征。因此，大速压条件下的静气动弹性预测应当加密计算工况点，且不宜用小速压条件下获得的变化规律外推至大速压范围。

4 结 论

本文基于RANS方程与静气动弹性平衡方程发展了一种CFD/CSD流固耦合静气动弹性计算方法。为了提升计算精准度、鲁棒性及效率，CFD采用了多块对接结构化网格进行分区并行计算，并利用多重网格技术来加快计算收敛速度，利用RBF结合TFI技术来生成结构化动网格，通过TPS方法来进行气动/结构的数据插值，并开发了相应代码给予实现。通过典型风洞试验结果验证了方法和代码的有效性。基于某大型飞机翼身组合体模型，数值模拟研究了跨声速时某典型大展弦比机翼的静气动弹性特性，获得了一些具有代表性的结论，可以提供给大型飞机气动/结构设计人员参考。

[1]Thomas G I. Unique testing capabilities of the NASA langley transonic dynamics tunnel, an exercise in aeroelastic scaling[R]. AIAA 2013-2625

[2]Mawilkinson D G, Blackerby W T, et al. Correlation of full-scale drag predictions with flight measurements on the C-141A aircraft-phase II, wind tunnel test, analysis, and prediction techniques[R]. NASA Technical Paper 2333, Feb. 1974

[3]Rodden W P, Johnson E H. MSC/Nastran Aeroelastic analysis user’s guide V68[M]. Los Angeles: MSC Software Corporation, 1994: 44-65

[4]Tian B Y. Computational aeroelastic analysis of airplane wings including geometry nonlinearity[D]. Graduate Education and Research University of Cincinnati, 2003, Paper for the degree of Doctor

[5]Zhang Z C, Liu F. Calculations of unsteady flow and flutter by an Euler and integral boundary-layer method on Cartesian grids[R]. AIAA 2004-5203

[6]方明祥, 白俊强等. 基于非结构动网格的跨音速颤振研究[J]. 机械设计与制造, 2007, 11: 11-19

[7]Livne E, Weisshaar T A. Aeroelasticity of nonconventional configurations: past and future[J]. Journal of Airplane, 2003, 40(6): 1047-1065

[8]Ramji Kamakoti, Wei S, et al. Time dependent RANS computation for an aeroelastic wing[R]. AIAA 2004-0886

[9]Robinson B A, Batina J T, et al. Aeroelastic analysis of wings using the Euler equation with a deforming mesh[J]. Journal of Airplane, 1991, 28(11): 781-788

[10]Guruswamy G P. Unsteady aerodynamic and aeroelastic calculations for wings using Euler equations[J]. AIAA Journal, 1990, 28(3): 461-469

[11]Rizzetta D P, Visbal M R, Gaitonde D V. Direct numerical and large-eddy simulation of supersonic flows by a high-order method[R]. AIAA 2000-2408

[12]Guarini S E, Moser R D, et al. Direct numerical simulation of a supersonic turbulent boundary layer at Mach 2.5[J]. Fluid Mech, 2000, 414: 1-33

[13]Weber C, Ducros F. A large eddy simulation of complex turbulent flows[R]. AIAA 98-2651

[14]李新亮, 傅德薰, 马延文. 可压缩均匀各向同性湍流直接数值模拟[J]. 中国科学A辑, 2002, 32(8): 716-724

[15]马汉东, 潘宏禄, 王强. 超声速平板边界层斜波失稳转捩过程研究[J]. 力学学报, 2007, 39(2): 153-157

[16]赵永辉. 气动弹性力学与控制[M]. 北京: 科学出版社, 2007

[17]Holger W. Computational aspects of radial basis function approximation[J]. Studies in Computational Mathematics, 2006, 2: 231-256

[18]Soni B K. Two- and three-dimensional grid generation for internal flow applications of computational fluid dynamics[C]//7th Computational Physics Conference, Cincinnati, OH, AIAA 1985-1526

[19]郦正能, 程小全, 方卫国. 飞行器结构学(第2版)[M]. 北京: 北京航空航天大学出版社, 2010.