大展弦比柔性机翼气动弹性分析中的气动力方法研究进展

杨 超, 杨 澜, 谢长川

(1. 北京航空航天大学 航空科学与工程学院, 北京 100191;2. 航空器先进设计技术工信部重点实验室, 北京 100191)

0 引 言

自20世纪90年代末起，由于长航时无人机、大型运输机和大型客机等长航时飞行器的高性能要求，高升阻比和轻质结构的设计充分体现在大展弦比机翼设计中，随之而来的一类新的非线性气动弹性问题开始受到关注，即大柔性飞行器的大变形几何非线性气动弹性问题。采用轻质材料的大展弦比机翼是该问题的主要研究对象，其力学本质在于结构求解中的小变形假设不再适用，结构受力变形后的平衡态相对未变形的结构呈现明显的几何差异，结构的承载和变形状态引起的几何非线性因素使得结构静、动特性发生改变，并且改变静、动气动弹性耦合关系，从而使气动弹性的研究及应用面临新的挑战。

几何非线性气动弹性的研究从理论方面与一般气动弹性力学主要有以下区别：其一为结构几何非线性理论,主要解决大变形情况下的结构静、动力学分析；其二为曲面气动力理论研究，主要解决结构大变形条件下，边界条件依赖于变形状态的定常和非定常气动力计算方法；其三为结构/气动界面耦合方法研究，主要研究适用于空间大变形的多维插值问题。

本文主要介绍大柔性机翼几何非线性气动弹性工程分析领域中气动建模方法的研究现状与进展，重点说明基于片条理论、面元法和CFD技术等气动建模方法的特点及其主要适用的问题，力图展示气动建模在几何非线性气动弹性研究中的应用特点，供相关气动弹性理论与工程应用、空气动力学等方面的研究人员参考。

1 片条理论

片条理论是气动弹性研究中使用最早且广泛应用的一种气动建模方法，其核心思想是利用二维流(无限翼展机翼)的简单结果来计算有限翼展升力面的定常/非定常气动力。将大展弦比升力面沿展向划分为若干窄条网格，基于二元翼段的气动理论计算每个片条的气动力，再将展向网格之间的气动干扰作用考虑在局部气动导数中进行三维修正，由此得到大展弦比机翼的气动特性。此类方法建立的气动模型简单，能与结构模型高效配合，并且能够考虑气动非线性作用，对于大展弦比机翼静、动气动弹性分析来说能够得到较为合理的计算结果，可用于飞行器气动弹性设计的初步阶段以及非线性机理研究。

基于片条理论计算二维翼型气动力的方法主要有Theodorsen方法、ONERA方法以及有限状态理论等，以下分别介绍。

1.1 Theodorsen方法

Theodorsen方法是一种基于线化理论的二维不可压流谐振荡非定常气动力理论，1934年由T. Theodorson提出[1]。在此基础上，建立了解析法精确求解二维翼面低速颤振问题的有效方法。Theodorsen理论能够得到简谐运动下精确的频域解析解，在20世纪早期的大展弦比飞机工程颤振计算中得到了广泛应用[2-4]。

谢长川等人将曲面修正的片条理论和Theodorsen理论相结合，建立了考虑机翼大变形的非定常气动力的计算方法，基于结构动力学准模态的思想进行大变形平衡态附近的气动弹性稳定性分析[5]。基于片条理论解决大变形曲面气动力问题时，重点是要考虑由于升力面几何形状引起的片条法向变化。气动力片条需在变形后的位置上定义，并建立片条局部坐标系。由结构分析得到的变形以及所有的气动力参数都从未变形机翼的总体坐标系XYZ变换到变形后的片条局部坐标系xyz中，如图1所示。

图1 变形前后片条坐标系示意图[5]Fig.1 Local strip coordinate systems before and after the deformation[5]

在变形后片条局部坐标系下，计算每个片条的定常迎角与非定常位移，再结合平衡位置附近的准模态信息，便可基于Theodorsen理论进行颤振分析。某大展弦比机翼的分析结果表明[5]，考虑几何非线性效应后，随着迎角增大，气动载荷增加，机翼变形增大，机翼颤振速度随之下降。如表1所示的基于曲面气动力的颤振速度结果可看出，机翼迎角为4°时，翼尖垂直位移达到半展长的15%，相对于无载荷状态，即不考虑几何非线性影响的计算结果，颤振速度下降了12.2%。由此可见，在大变形机翼的颤振分析中，几何非线性因素不可忽略。对比表1所示的两组不同气动力方法得到的颤振速度结果可以看出，两种方法采用同样的结构模态数据，但分别使用变形后的曲面气动力和未变形机翼平面气动力，在4°迎角时两组颤振结果偏差较大，说明大变形曲面气动力效应对颤振结果有较大的影响。尤其是在变形较为显著时，颤振分析应当采用曲面非定常气动力方法进行计算。

表1 考虑几何非线性的机翼颤振特性[5]Table 1 Flutter results considering geometric nonlinearity[5]

1.2 ONERA方法

20世纪80年代法国宇航研究院的Tran与Petot等人为研究柔性直升机旋翼在大迎角非定常动力失速情况下的气动弹性特性，首次提出ONERA模型[6]，它是一个半经验、非定常、非线性的二维气动力计算方法。ONERA模型中的方程系数由风洞试验数据拟合得到，其具体取值依赖于具体的试验状态和飞行条件，能够较准确地预测翼型的气动特性。

为考虑动失速的影响，向锦武等人基于ONERA动失速气动力模型与片条假设，结合几何精确完全本征运动梁模型，建立了大展弦比柔性飞机非线性气动弹性与飞行动力学耦合模型[7]。如图2所示，大展弦比机翼沿展向划分为若干互相独立的片条气动单元，通过考虑翼型相对梁参考线的偏移与扭转，应用ONERA气动理论得到全机的气动载荷。其研究结果图3表明，当机翼变形较大时，失速首先发生于翼尖而后向翼根扩展，失速范围有限且全机升力损失相对较小。

图2 柔性飞机非线性气动弹性与飞行动力学耦合建模示意图[7]Fig.2 Coupled nonlinear aeroelastic and flight dynamic modelling of flexible aircraft[7]

图3 飞翼布局柔性飞机满载时局部有效迎角的时间历程[7]Fig.3 Time history of local effective α of flexible flying wing in full load case[7]

赵永辉等人采用简化梁和ONERA模型研究了大展弦比机翼的大迎角颤振特性[8]；张亮等人采用Dowell-Hodges梁和ONERA模型，并使用线性化方法分析了大展弦比机翼的颤振特性[9]。王伟等人基于共旋有限元(Co-rotational，CR)理论和ONERA模型，提出了一种适用于大柔性机翼的非线性气动弹性求解方法[10]。其研究结果表明，非定常气动力的动态失速效应限制了非线性气动弹性响应振幅的增加，但也会带来疲劳问题，结构设计时需要特别注意。

1.3 有限状态理论

有限状态理论于20世纪90年代由Peters建立[11-12]，针对不可压、二维薄翼可直接由势流方程推导出，对翼型的运动状态没有限制，气动力状态以诱导流的系数形式给出。有限状态气动(入流)理论可以只考虑一定数目的状态和频率范围，无需中间的数值分析。一阶状态方程很容易与控制方程耦合，可用于频域、Laplace域和时域分析。与Theodorsen、Wagner函数相比，有限状态方法只需要很少的状态就能得到很好的精度。

Su等人采用应变基几何非线性有限元方法与有限状态理论建立了一种非线性气动弹性分析方法[13-14]。有限状态(入流)模型基于薄翼假设，能够计算大变形或大幅运动下的二维非定常气动力。在与应变基有限元方法耦合时，局部气动片条坐标系需要根据瞬态的结构节点速度与位移向总体坐标系转换，从而描述精确的几何信息。

Palacios等人研究了不同的结构、气动建模方法在大柔性飞行器的非线性飞行动力学中的应用[15]。其中，涡格法和片条理论都能够用于几何非线性气动弹性问题中的非定常气动力建模。图4给出了基于有限状态方法和Theodorsen理论得到的二维非定常气动力的对比，可以看出，在一定的减缩频率范围内，有限状态方法能够与Theodorsen理论达到相近的精度，但是当选取的状态过多时，由于有限状态理论中对方程的近似处理[16]，其计算精度会下降。其研究还对比了基于片条理论(分别使用了有限状态方法Finite-state，和指示函数法Indicial Response)和非定常涡格法(Unsteady Vortex-Lattice Method，简称UVLM)计算的三维非定常气动力，如图5所示。结果表明，在较低的减缩频率下，翼尖的三维效应对非定常气动力的影响很大，尤其对小展弦比机翼，片条理论需要进行翼尖修正才能得到较好的结果。

此外，Cesnik等人也基于有限状态理论建立了非线性气动弹性仿真工具箱——UM/NAST，并进行了一系列关于高空长航时(High Altitude Long Endurance，简称HALE)飞机的几何非线性气动弹性研究，包括非线性气动弹性建模[17]、机动载荷研究[18]、颤振边界扩宽[19]以及阵风载荷减缓研究[20]等。

(a) 幅值

(b) 相位

(a) AR=2, k=0.05

(b) AR=10, k=0.05

(c) AR=2, k=0.5

(b) AR=10, k=0.2

2 面元法

面元法作为一种气动力计算的工程方法，近40年来一直发展不衰。面元法将复杂组合体三维绕流问题简化为一个二维积分方程，从而大大降低数值计算量[21]。现今气动弹性工程分析中的面元法种类众多，其中，三维偶极子格网法(Doublet-Lattice Method，DLM)是气弹领域中最常用的一种频域气动力计算的工程方法[22]。涡格法(Vortex-Lattice Method，VLM)在时域非定常气动力的计算方面有突出的优势。

传统的气动弹性分析中的面元法不考虑机翼升力面的曲面效应，以MSC.Nastran为代表的商业软件基于平板气动力进行分析和计算，对于传统线性问题具有足够的精度。然而对于几何非线性问题，由于大柔性机翼在外载荷的作用下发生较大的弹性变形，机翼升力面形状与原本的未变形状态产生显著差别，曲面气动力与结构的相互作用使得曲面气动力计算成为不可避免的气动弹性分析部分。

对于大柔性飞机的曲面气动力计算，面元法中边界条件与传统线性边界条件有显著区别，其核心问题是考虑气动构型变化的几何精确的曲面边界条件，边界条件方程可简单表达为[23]

(1)

其中v为物面运动速度，n为物面法向。

本文主要对曲面偶极子格网法和定常/非定常曲面涡格法的研究进展进行介绍。表2总结了三种曲面气动力方法的特点与适用范围。其中曲面涡格法能够较好地对三维曲面气动面进行建模，并且模型简单、计算效率高，在大柔性飞机曲面时域气动力计算方面有明显优势。曲面偶极子格网法延续了平面偶极子格网法的求解思路，但是曲面气动面的建立使其与大柔性飞机气动力计算需求相一致。作为频域气动力计算方法，曲面偶极子格网法易与颤振方程结合，便于气动弹性颤振问题的求解。

表2 三种曲面气动力方法对比Table 2 Comparison of three nonplanar aerodynamic methods

2.1 曲面偶极子格网法

偶极子格网法(DLM)是目前气动弹性工程分析中最为常用的线性非定常气动力计算方法，目前商用软件(如MSC. Nastran)中的非定常气动力模块就采用了偶极子格网法[24]。其基于小扰动线化位势流频域方程，在机翼表面布置压力偶极子形式的基本解，避免了对尾流区的处理，是工程颤振分析上较常用的一种非定常气动力计算方法[25]。

传统基于气动力平面假设的偶极子格网法模型不能考虑机翼的弯曲及扭转变形的影响，对于结构发生大变形的柔性机翼是不准确的。为建立适用于大变形机翼的非定常气动力计算，谢长川等人对其进行三维曲面修正，其核心问题是修正曲面的边界条件，使用更广泛的物面边界条件[23]：

(2)

图6 曲面网格偶极子线及法洗方向[27]Fig.6 Nonplanar doublet lattice and normal wash[27]

Patil和Hodges将非线性梁模型与曲面偶极子格网法结合做了简单的应用研究[28]。从工程分析的实际需要来看，曲面偶极子格网法能够有效地考虑大弯曲和小扭转变形情况下升力面的频域非定常气动力，针对可以引入动态小扰动假设条件的机翼部件、全机气动弹性稳定性与响应问题可以给出很好的工程近似，便于继承现有的传统线性气动弹性分析工具。

谢长川等人基于曲面偶极子格网法开展了几何非线性气动弹性稳定性分析[27]。如图7所示，由于不同迎角下的机翼静平衡构型的差异，考虑几何非线性效应的颤振分析明显比传统线性结果更接近风洞试验数据。对于大变形情况下的机翼气动弹性分析，由于颤振速度与静平衡状态的变形与应力状态密切相关，因此基于静变形平衡构型建立气动模型的曲面偶极子格网法是一种十分有效的气动建模方法。

图7 不同迎角下的颤振速度[27]Fig.7 Flutter speeds under different angle of attack[27]

2.2 定常/非定常曲面涡格法

对于薄机翼，涡格法(VLM)也是一种十分有效的气动力计算方法。涡格法由势流方程推导而来，将机翼离散成沿弦向和展向分布的附着涡并在后缘布置有自由涡线来模拟尾流区。曲面涡格法能考虑机翼弯度的影响，且计算量小，并保证一定的精度，在气动弹性问题中也具有突出的优势。

涡格法在定常、非定常气动力分析方面都十分有效。图8和图9分别展示了定常、非定常曲面涡格法的网格布置情况[23]。定常与非定常曲面涡格法求解思路一致，但在具体尾涡流场建模和气动力计算方面有所差别。在定常涡格法中，翼面自由涡由后缘涡格拖出，平行于来流方向。而非定常涡格法的尾流区仍利用涡环进行模拟，随时间步推进。在不同时间计算步中，机翼表面涡环环量不同，尾流区的涡环分布也随之变化，呈现出显著的非定常特点。

图8 定常曲面涡格中涡环单元布置情况[23]Fig.8 Vortex rings in nonplanar steady VLM[23]

图9 非定常曲面涡格法气动模型[23]Fig.9 Aerodynamic model of nonplanar UVLM[23]

刘燚、谢长川等基于定常曲面涡格法进行快速的静气动弹性分析[29-30]。非定常曲面涡格法(UVLM)作为一种时域气动力计算方法，由于建模简单、计算效率高、易与结构求解结合，成为流固耦合求解常用的工程气动力计算方法[31]。UVLM能给出机翼尾流区的位置并且气动节点随结构节点不断变形，适应了柔性机翼非线性曲面气动力计算的需求，因此在柔性飞机上得到了广泛的应用。

P. C. Chen用非定常涡格法结合考虑弯曲和扭转的几何非线性有限元梁单元(Hodges内蕴梁)，在时域范围内对控制方程进行积分，讨论了HALE机翼的动力学响应问题[32]。他还基于UVLM研究了离散阵风气流分离、机翼各种迎角范围内的气动力模型并直接进行时域的动力学响应分析[33]。Cesnik将UVLM与非线性板单元结合起来，加入了失速模型，并研究了各向异性材料对平飞和颤振的影响[34]。刘燚等基于非定常曲面涡格法研究了大展弦比机翼的阵风响应[23]。如图10所示，非定常涡格法的自由尾涡模型能够模拟机翼尾涡的非定常规律。图11为翼尖垂直位移响应的试验数据与计算结果对比，可以看出基于UVLM的气动弹性响应时域仿真结果与风洞试验数据吻合得很好。

图10 非定常涡格法自由尾涡模型[23]Fig.10 Free wake model in UVLM[23]

Ritter等对比了基于片条理论和非定常涡格法等不同气动力方法对大展弦比柔性飞行器的非线性静、动气动弹性分析[35]。如图12所示，文献[35]分别基于三种气动弹性求解器(UM/NAST，DLR和NASTRAN)计算了全机的阵风响应问题。由于片条理论不考虑不同气动面之间的相互诱导作用，基于片条理论(图12中UM/NAST strip theory)与基于非定常涡格法(图12中UM/NAST UVL, DLR UVL)的气动弹性响应计算结果存在较大偏差。

图11 翼尖垂向位移响应的试验数据与计算结果对比[23]Fig.11 Analysis and test results ofvertical wing tip deflections (F=2 Hz)[23]

(a) Wgust=0.2 m, right wingtip

(b) Wgust=1.0 m, right wingtip

3 CFD技术

随着计算机水平的大幅提高，CFD技术已经被广泛应用于各种气动弹性机理研究和工程计算中[36]。近年来的气动弹性计算中也逐渐采用了以跨声速小扰动方程、Euler方程或N-S方程为基础的CFD技术来计算非定常气动力[37-39]。这种方法直接从流动的基本方程出发，使用的假设相对较少，模拟了流动的本质特性，可以反映出气动力的非线性特性[40-42]。但基于CFD技术的气动弹性分析计算量过大，工程上一般仅用于典型或严重工况下的校核，而在设计初期时，仍普遍采用基于面元法的快速分析方法进行大量工况的计算。

杨国伟对近年来计算气动弹性力学研究方面的若干进展进行了综述[36]，并指出计算气动弹性力学的研究重点在于两方面：一是计算流体力学(CFD)和计算结构动力学(Computational Structural Dynamics，CSD)及其耦合计算方法的研究；另一方面是发展通用性好、计算效率高、适用性强的动网格技术和流体、结构间的精确数据插值方法。对于大变形柔性机翼来说，基于CFD技术的气动弹性分析，其研究重点也同样在于这两方面的突破。尤其是结构大变形对动网格的技术要求更高，如何在显著大变形情况下生成正交性和光顺性良好的网格成为研究难点之一。

徐敏等人针对大展弦比柔性机翼，发展了一种基于CFD/CSD弱耦合的非线性气动弹性计算技术，研究了非线性静气动弹性的分析方法[43]。崔鹏等人将非线性有限元和CFD计算相结合，采用更新的Lagrange方法分析了结构大变形引起的几何非线性气动弹性问题，并研究了切尖三角翼的极限环(Limit Cycle Oscillation, LCO)现象[44]。Patil 和Hodges将非线性气动力(CFD)与精确梁理论以松耦合的方式研究了大展弦比机翼的静气弹行为，着重研究了非线性气动力对弹性效用的影响及其重要性，研究发现CFD计算气动力比面元法偏小，因此低阶面元法对发散和颤振的计算可能过于保守[45]。Palacios等基于欧拉方程和几何非线性梁模型，开展了紧耦合下的气动弹性响应研究，分析了大展弦比飞机的稳态飞行问题[46]。在此基础上，Hallissy等开发了一种针对大柔性飞行器的多学科高精度气动弹性仿真工具[47]，其研究结果表明，在静气动弹性分析方面，基于CFD/CSD耦合的高精度方法与基于片条理论的低精度方法的分析结果吻合得很好，如图13所示。但在颤振预测方面，基于低精度方法的颤振速度预测对算例机翼较为保守，比基于CFD的高精度方法得到的颤振速度高约2%～5%。

虽然近年来基于CFD/CSD耦合方法的气动弹性研究取得了较大进展，能够解决机翼部件的非线性静、动气动弹性问题，但是计算规模大，效率低，难以满足工程气动弹性设计的需求。特别是对于大柔性飞行器的全机气动弹性问题，基于CFD技术的气动力分析在与飞行力学、结构动力学及飞行控制等的多学科耦合方面仍存在一定困难。

为提高计算效率，基于CFD技术的非定常气动力降阶模型(Reduced-order Model，ROM)是一种十分有效的传统对象的气动弹性计算方法方法。然而现有的气动力降阶模型大多不能考虑机翼大变形情况下的气动力降阶问题；对于大变形机翼的几何非线性气动弹性问题，几乎没有基于ROM的相关文献研究。

图13 静气弹位移结果对比[47]Fig.13 Comparison of static aeroelastic displacements[47]

4 结论及建议

大柔性飞行器由于其结构柔性大，重量轻，在飞行过程中易发生大变形而引起几何非线性问题。本文回顾和总结了大柔性机翼几何非线性气动弹性研究和工程领域主要使用的气动建模方法及其应用情况，着重介绍了基于片条理论、面元法和CFD技术的非线性气动弹性力学的研究进展情况。在几何非线性气动弹性领域中，气动建模方法的研究重点在于解决结构大变形条件下边界条件依赖于变形状态的定常和非定常气动力计算方法，与常规线性气动弹性分析显著不同。通过分析国内外在几何非线性气动弹性领域的研究进展，针对此类问题主要有以下几方面的结论和建议：

(1) 对于静气弹问题，基于片条理论或面元法的几何非线性气动弹性分析基本能够满足精度要求，并且建模简单、计算效率高；基于CFD/CSD耦合的静气弹分析能够保证更高的计算精度，但计算规模相对较大，建模较为复杂。因此，对于初步设计阶段的大量工况分析可基于片条理论或面元法进行快速气弹分析；对于典型或严重工况，可以采用CFD/CSD耦合方法进行高精度分析。

(2) 对于动气动弹性稳定性，即颤振问题的研究，大变形情况下的非线性颤振分析依赖于给定工况下的静平衡构型的应力和变形信息，非线性颤振速度与线性分析结果有显著差异。因此，对于大变形问题，应当重视几何非线性对结构频率和气动构型的影响，采用非线性颤振求解方法进行分析。对于大变形机翼部件的颤振计算，基于片条理论、面元法和CFD技术的非线性颤振求解方法都能够得到较为理想的结果。但对于大变形机翼全机非线性颤振问题，基于CFD/CSD耦合的求解方法目前仍存在一定的困难，需要解决与飞行力学、结构动力学及飞行控制等的多学科耦合计算方法和效率问题，才能在工程中应用。

(3) 对于气动弹性动响应问题，基于片条理论的气动弹性分析方法往往不够准确。当迎角较小、不考虑失速问题时，非定常曲面涡格法能够十分有效地进行大变形情况下的时域响应计算。但对于较为复杂的流场，如气动力具有跨声速、失速等明显的非线性特性时，应当进行CFD/CSD耦合时域仿真。此外，对于大变形机翼全机的气动弹性响应问题，基于CFD/CSD耦合的求解方法目前仍在许多方面存在困难，仍然需要解决与飞行力学、结构动力学及飞行控制等的多学科耦合计算方法和效率问题。

(4) 目前基于CFD技术的非定常气动力降阶模型在大变形几何非线性气动弹性研究方面的进展很少，亟需发展适用于大展弦比柔性机翼部件和全机气动弹性分析的高效高精度的基于CFD技术的气动力降阶分析方法。

为适应高空长航时无人机、大型运输机和大型客机等长航时飞行器的研制发展需求，适用于大柔性飞机气动弹性分析的气动力建模方法仍在很多方面需要发展和完善。更高效、更高精度的气动力计算方法始终是气动弹性研究领域的挑战，其气动建模技术的发展也必将推动气动弹性分析水平的不断提高。