基于高阶谐波平衡的跨声速颤振高效预测方法

刘南，郭承鹏，白俊强

1. 中国航空工业空气动力研究院 气动研究与试验二部，沈阳 110034 2. 高速高雷诺数气动力航空科技重点实验室，沈阳 110034 3. 西北工业大学 航空学院，西安 710072

颤振是一种典型的气动弹性不稳定性现象，可能导致飞行器的结构遭到灾难性破坏，在飞行包线内绝对不允许发生颤振现象且留有一定的裕度。因此，准确预测飞行器颤振特性对结构设计和强度校核具有重要意义，传统工程颤振计算采用基于谐振荡假设的偶极子格网法，在跨声速阶段该方法的精度大幅下降。然而，飞行器的跨声速颤振速压边界远低于其他速域，而且激波运动及附面层分离常引发各种复杂的非线性颤振现象，显著地影响飞行安全和任务性能。

针对飞行器跨声速颤振问题，计算流体力学(Computational Fluid Dynamics，CFD)方法能够较准确地预测跨声速流场，因此耦合非定常CFD方法和计算结构动力学(Computational Structural Dynamics，CSD)成为跨声速颤振分析的重要手段[1]。但是，三维流场的空间自由度通常达到106量级以上，利用CFD方法计算非定常气动力需要耗费大量的计算资源和计算时间。例如波音公司的Hong等[2]利用CFD/CSD耦合方法分析HSCT构型的颤振特性时，获得每个状态的时域响应需要花费长达108个CPU小时。

为了提高分析效率，研究人员提出了若干种气动弹性高效分析手段，主要有Volterra级数、本征正交分解和谐波平衡(Harmonic Balance, HB)方法[3]，其中HB方法因其物理意义明确且适用于非线性问题等优点受到广泛关注[4]。Hall等[5]进一步改进HB方法，不直接求解流场解的谐波系数，而是通过离散Fourier变换和逆变换建立谐波系数和一个周期内等距分布的流场变量之间的关系，对后者进行求解，不仅提高了分析效率，而且实现较为简单，这种方法通常被称作高阶谐波平衡(High-Order Harmonic Balance, HOHB)方法。该方法公式中显式地出现运动频率，因此对于频率已知问题，HOHB方法表现优异。Ekici等[6]利用无黏Euler求解器Duke ROTOR建立HOHB分析方法，从而对旋翼前飞和悬停状态进行分析，表现出较高的可靠性。Woodgate和Badcock[7]构造了一种隐式HOHB求解器，对NACA0012翼型和F-5机翼周期性强迫俯仰运动进行分析，与非定常计算相比，在保证精度的前提下，所需CPU时间降低10倍以上。陈琦等[8]对现有计算程序进行改造，实现HOHB方法，对跨声速翼型和超声速钝锥周期性强迫俯仰振荡进行分析，结果表明，二阶HOHB方法能够得到时域响应与非定常计算非常接近，计算时间降低80%。Ronch等[9]采用HOHB方法分析了跨声速DLR-F12翼身组合体的动导数，并与线性频域法进行对比，HOHB表现出较高的准确性。贡伊明等[10]采用带预处理的广义极小残差算法提高了HOHB方法的计算收敛性。

对于运动频率未知问题，例如圆柱绕流的漩涡脱落频率、颤振频率、极限环振荡频率等，需要添加额外的条件。其中一种是相位固定条件，通过相位对频率的梯度进行分析。利用这种方法，Thomas等[11-12]基于HOHB方法建立了一种非线性极限环求解器，对无黏NACA64A010翼型和黏性NLR7301翼型跨声速极限环振荡开展研究。Liu和Dowell[13]研究了HDHB方法在带操纵面极限环非线性的二维翼型极限环振荡问题中的应用。刘南等[14]研究了高阶谐波平衡方法中的非物理解来源，通过扩充非线性项的子时间层消除了非物理解且降低了计算所需的阶数，并将之应用于二维两自由度立方非线性极限环振荡问题中。但是相位固定条件对于复杂问题则非常繁琐和耗时。Ekici和Hall[15]提出了一种基于残差L2范数最小的频率更新方法，在单自由度非线性极限环振荡的数值模拟中得到较好的结果。Yao和Marques[16]指出Ekici方法在多自由度问题中收敛性变差，并提出改进措施，用于分析二维翼型和三维机翼的无黏极限环振荡。但是基于残差L2范数最小的频率更新方法较为依赖给定的频率初值，可靠性较差。

针对颤振频率未知的问题，利用HOHB方法结合结构振型首先快速地获得广义气动力影响系数(generalized Aerodynamic force Influence Coefficient, AIC)矩阵，从而获得系统阻尼和来流速度之间的关系，建立高效颤振分析流程，避免了在HOHB方法中进行频率求解的难题，有效地提高了颤振分析效率。

1 颤振时间推进方法

基于CFD/CSD耦合的颤振时间推进方法涉及到非定常流场求解、结构有限元求解以及气动/结构数据传递等。介绍本文所采用的流场求解、结构运动方程求解、数据传递和流程建立等部分，最终利用AGARD445.6国际标摸验证分析方法的精度和可靠性。

1.1 CFD求解方法

计算坐标系下无源项的三维可压缩RANS方程可以写成：

(1)

颤振计算涉及到气动模型的几何变形，在非定常计算物理时间步后均需调用网格变形子程序更新空间网格。目前，研究人员提出多种网格变形方法，但是这些方法在变形效率、变形能力和变形后网格质量方面表现出明显的优缺点[21]。

因此，本文利用背景网格的思想，结合径向基函数[22](Radial Basis Function，RBF)和无限插值[23](Trans Finite Interpolation，TFI)两种网格变形方法的优点，建立了一套高效高鲁棒性网格变形方法[24]。

1.2 结构有限元求解

结构运动方程可以统一写成：

(2)

式中：M、D和K分别为结构质量阵、阻尼阵和刚度阵；q为结构变形；f为气动力。采用Patran对结构进行建模，通过Nastran SOL103求解器可以获得结构振型φm,m=1,2,…及对应的特征频率ωm,m=1,2,…。仅考虑前M阶结构振型，组装成矩阵：

(3)

结构变形q可以通过结构振型与广义位移的线性叠加得到，

(4)

利用结构振型可以建立广义结构运动方程：

(5)

式中：广义质量阵、阻尼阵和刚度阵的表达式为

(6)

广义力的表达式为

(7)

1.3 气动/结构数据传递

能量守恒是气动/结构交界面数据传递需要满足的一种重要条件[22]，利用虚功原理得到

(8)

式中：δW为虚功；δu为虚位移；f为力向量，下标a代表气动节点，s为结构节点。定义一个耦合矩阵，气动节点和结构节点的位移之间的关系为

ua=Hus

(9)

通过式(8)和式(9)，力插值公式可以写成：

fs=HTfa

(10)

通过RBF方法[22]可以构造矩阵H，还可以满足力和力矩守恒性、插值外形光滑性等要求，此处不再赘述。

对于结构振型，将结构振型利用矩阵H插值到气动表面得

Φa=HΦs

(11)

式中：Φa、Φs分别是结构振型在气动物面网格点和结构网格点上的表达。式(11)结合式(4)和式(7) 即可满足能量守恒要求。因此，仅需将结构振型插值到气动表面网格上获得Φa作为计算输入文件，在各个物理时间步通过式(12)和式(13) 获得气动力和位移：

(12)

(13)

1.4 时域计算流程

将广义结构运动方程转化为状态空间方程形式：

(14)

式中：

(15)

图1 气动弹性时域分析流程Fig.1 Process of aeroelastic time marching analysis

颤振时间推进方法的流程如图1所示，图中展示了流场和结构的时间推进和数据传递路线，是一种典型的串行处理方法。

1.5 标模验证

AGARD445.6机翼于1961年在NASA兰利中心TDT(Transonic Dynamics Tunnel)风洞进行试验[25]，其中三号软模型成为颤振分析程序验证的标准模型[26-27]。利用本文建立的颤振分析方法预测AGARD445.6机翼的无黏和黏性颤振边界，试验结果与本文计算结果的对比如图2所示，其中EXP、Inviscid和Viscous分别代表试验结果、无黏计算结果和黏性计算结果，横坐标为来流马赫数，纵坐标Qinf和Freq分别代表来流动压和振动频率。由图可见：本文所建立颤振时间推进方法得到的黏性结果与试验结果吻合较好，无黏结果在低超声速工况颤振速压偏高，与参考文献[26-27]十分吻合，表明了本文方法的可靠性。

图2 AGARD445.6机翼颤振边界结果对比Fig.2 Comparison of flutter boundary result of AGARD445.6 wing

2 高效颤振分析方法

2.1 HOHB流场求解方法

式(1)可以写成

(16)

其中

(17)

如果已知式(16)的解具有周期性且基频为ω，则Q(t)可以通过Fourier级数展开：

(18)

式中：NH为谐波阶数。将式(18)对时间求导：

(19)

式(16)方程右端残差项R的各阶Fourier系数可以通过式(20)～式(22)推导或数值计算获得：

(20)

(21)

(22)

再根据方程式(16)左右两侧相等的原则，可以得到

(23)

写成矩阵形式：

(24)

其中

(25)

其中

Q(t2NH)]T

为等距分布的各子时间层ti,i=0,1,2,…,2NH上的时域解，ti分布为

(26)

矩阵E为离散Fourier变换矩阵：

(27)

将式(25)代入式(24)，得到：

(28)

式(28)两边同乘以E-1：

(29)

直接将式(29)中第2项简化为

(30)

其中

R(t2NH)]T

将式(30)代入式(29)，可得

(31)

式中：U=E-1JE。引入伪时间导数项驱动式(31) 残差趋于0，构造求解方程：

(32)

2.2 频域结构运动方程

对于颤振等周期性问题，模态位移可以表示为

(33)

(34)

将式(33)和式(34)代入结构运动方程式(5)可得：

(35)

(36)

其中

假设结构模态位移的运动形式为

(37)

式中：p=g+jω，ω为振荡频率(为了简单起见，本文没有采用减缩频率k，而是直接采用频率ω)，g为系统阻尼，是一种更准确的物理衰减率的近似。将式(37)代入广义结构运动方程式(5)可得

(38)

根据式(38)，可得

(39)

(40)

其中式(39)是系统静变形方程，可以获得系统静平衡位置，式(40)就是颤振频域控制方程。将式(36)代入(40)，可得

(41)

参考Chen等提出的g方法，引入变量：

(42)

将上述变换代入式(41)，可得

(43)

综合式(42)和式(43)得到

Ly=gy

(44)

式中：

L=

至此，就通过广义力影响系数矩阵将颤振边界问题转化为复特征值求解问题。具体的分析流程如下：

1) 利用HOHB求解器快速获得一系列频率{ωi,i=1,2,…,Nω}下的广义气动力影响系数矩阵F(jω)。

2) 在某频率ωi下，得到一系列来流速压Qi,i=1,2,…,NQ下的矩阵L，并采用复矩阵特征值分析得到不同来流速压下的特征值。

3) 根据特征值随来流速度的变化趋势，得到特征值虚部为0时特征值实部g和来流速压Q。

4) 得到特征值实部g和来流速压Q随频率的变化趋势，当特征值实部g=0时，对应的ω和Q分别就是颤振临界频率和动压。

上述就是本文建立的跨声速颤振高效预测方法。综上可见，本文避免了直接求解颤振频率ω的难题，而是利用HOHB方法获得多个频率下的广义气动力系数矩阵，通过频域法计算颤振边界。

3 二维翼型颤振边界预测

对于仅具有沉浮和俯仰两自由度的二维翼型，广义力和模态位移之间的关系为

(45)

图3 Ma=0.80时通过强迫简谐俯仰运动得到的广义气动力影响系数Fig.3 Generalized aerodynamic force influence coefficients obtained by harmonic pitch motion at Ma=0.80

图4 Ma=0.85时通过强迫简谐俯仰运动得到的广义气动力影响系数Fig.4 Generalized aerodynamic force influence coefficients obtained by harmonic pitch motion at Ma=0.85

Ma=0.80和0.85时沉浮和俯仰模态阻尼系数随来流速度的变化趋势如图5所示，横坐标是来流速度，纵坐标是模态阻尼系数。在Ma=0.80时，沉浮模态是颤振临界模态，黏性对阻尼曲线影响较小；在Ma=0.85时，无黏结果在V=33.5 m/s附近沉浮模态的阻尼系数变成正数，系统发生颤振，然后在V=140 m/s附近阻尼系数又重新变成负数，系统又回到稳定状态，最后在V=148.8 m/s 附近俯仰模态的阻尼变成正数，系统再次发生颤振，共有3个临界点；而黏性结果仅在V=128 m/s附近沉浮模态阻尼变成正数，表明系统发生颤振。

本文建立的HOHB高效颤振分析方法与时间推进方法得到的无黏和黏性颤振边界对比分别如图6和图7所示，其中Time marching和HOHB分别是时间推进和HOHB两种方法的计算结果，横坐标是来流马赫数，图6(a)和图6(b)中纵坐标分别是颤振速度和频率。由图可见，本文方法具有较高的预测精度，能够准确捕捉无黏的“S”形颤振边界[30]、高马赫数(Ma≥0.90)时的“锁定”现象[31]以及黏性结果中Ma≥0.85时由激波诱导附面层分离导致的颤振速度增大现象[32]。

图5 Ma=0.80和0.85时模态阻尼系数随来流速度的变化趋势Fig.5 Variation of modal damping coefficients with freestream velocity at Ma=0.80 and Ma=0.85

两个方法计算所用时间对比如表1所示。由表可见，HOHB方法所用时间约为时间推进方法的1/6，但是值得注意的是表中估算是较为保守的，仅采用5次非定常计算很难得到较为准确的颤振边界，而计算10个频率下的广义力影响系数矩阵偏多。

图6 时间推进方法和HOHB法计算得到的无黏颤振速度和频率对比Fig.6 Comparison of inviscid flutter velocities and frequencies obtained with time marching and HOHB methods

图7 时间推进方法和HOHB法计算得到的黏性颤振速度和频率对比Fig.7 Comparison of viscous flutter velocities and frequencies obtained with time marching and HOHB methods

方法计算时间加速比时间推进5×33.95=169.75HOHB10×2×1.43+0.58=29.185.82

4 三维机翼颤振边界预测

然后以三维Goland+机翼为研究对象，具体气动外形和结构参数参见Beran等[33]。选取前四阶结构模态(包括一阶弯曲模态、一阶扭转模态、二阶弯曲模态和二阶扭转模态)，频域内广义力和模态位移之间的关系为

(46)

令各个结构模态做一系列频率的简谐运动qm(t)=0.001sin(ωit),m=1,2,3,4，从而获得广义力影响系数矩阵，然后利用复特征值分析得到模态阻尼随来流速度的变化趋势，如图8～图10所示，其中Mode1、Mode2、Mode3和Mode4分别是前1至4阶结构模态阻尼系数，横坐标是来流速度，纵坐标是模态阻尼，分别是无外挂物构型的无黏阻尼系数矩阵以及外挂物构型的无黏和黏性阻尼系数矩阵。结果表明：Ma=0.80时颤振临界模态是一阶弯曲模态，而Ma=0.92时颤振临界模态由一阶弯曲模态转移至一阶扭转模态，这是由于激波引起压心后移导致的。

图8 无外挂物Goland+机翼无黏阻尼系数随来流速度的变化趋势Fig.8 Variation of inviscid damping coefficients with freestream velocity of Goland+ wing without store

图11至图13是颤振速度和频率对比，其中Time marching是时间推进方法计算结果，HOHB nmodes=4和HOHB nmodes=2分别是考虑前四阶模态和仅考虑前两阶模态的HOHB方法计算结果，横坐标是来流马赫数。在时间推进方法中，结构模态数量对计算效率基本没有影响(在颤振时间推进流程中广义结构运动方程的计算消耗相比非定常流场计算微乎其微)，但是在HOHB方法中需要使各个结构模态做简谐运动，

图11 无外挂物Goland+机翼无黏颤振边界对比Fig.11 Comparison of inviscid flutter boundaries of Goland+ wing without store

图12 带外挂物Goland+机翼无黏颤振边界对比Fig.12 Comparison of inviscid flutter boundaries of Goland+ wing with store

图13 带外挂物Goland+机翼黏性颤振边界对比Fig.13 Comparison of viscous flutter boundaries of Goland+ wing with store

从而获得广义力影响系数矩阵，计算效率与模态数成反比。结果表明：① 考虑所有结构模态的HOHB法结果与时间推进方法吻合很好，仅考虑前两阶模态也能够较好地预测跨声速颤振边界的变化趋势；② 对于不同的构型，高阶模态的影响有所不同，例如无外挂物构型，忽略二阶弯曲和扭转模态后HOHB法预测的颤振速度平均误差在1.0%以内，而带外挂物构型平均误差在3.0%左右；③ HOHB法所用时间仅为时间推进方法的1/6，能够显著提高颤振边界预测效率。

5 结 论

1) 基于高精度CFD方法和结构模态运动方程耦合求解，建立了颤振时间推进方法，通过对AGARD445.6机翼标摸进行测试，计算与试验结果吻合较好，表现出较高的计算精度和可靠性。

2) 对于颤振等频率未知的问题，很难利用HOHB方法直接求解。采取一种间接的方式——利用HOHB方法快速地获得广义气动力系数矩阵，进而求解颤振边界，建立高效颤振频域分析流程。

3) 基于HOHB建立高效颤振频域计算方法，并将之应用于二维和三维跨声速颤振边界预测，与时间推进方法计算结果吻合较好，表现出较高的效率和可靠性。