基于降阶模型的不同厚度翼型颤振边界预测

周 萌，高国柱

(中国电子科技集团公司第三十八研究所 浮空平台部，合肥 230088)

0 引言

为了解决CFD/CSD计算效率低下的问题，在20世纪90年代，杜克大学的Dowell[1]和美国国家航天航空局的Silva[2]等人将降阶模型(ROM)方法应用到气动弹性领域中。构造降阶模型的目的主要有两点[3]：一是降阶模型技术大大减小原数值模型系统的阶数，且在保证较为精确的系统动力学特征描述的同时大大减少计算耗时；二是降阶模型为学者们研究系统动力学特征提供了强有力的工具。目前，降阶模型主要分为两类：一类是对流场特征结构进行分析的降阶模型，其中以基于本征正交分解(Proper Orthogonal Decomposition，POD)技术[4]的降阶模型和基于谐波平衡(Harmonic Balance，HB)方法[5]的降阶模型为代表；另一类是基于系统外部输入和输出响应的系统辨识方法的降阶模型，这种方法以Silver的Volterra模型[6]、基于代理模型[7]和基于时间序列的非线性带外部输入的自回归滑动平均模型[8](NARMAX)的降阶模型为代表。HB方法是一种在频域内求解非线性周期流场的方法。POD方法是一种获得高维空间在低维空间内近似表达的方法，通过CFD计算快照构造相关矩阵，获得特征模态形成降阶子空间，将控制方程投影到系统特征来建立ROM。HB方法和POD方法是在全阶系统的基础上，建立新的控制方程并对原系统数值求解源代码进行修改，进而提高分析效率。系统辨识方法相较于HB方法和POD方法处理方式简单，适用性强，能够适用于线性和非线性问题。

目前，国内外研究人员基于不同的ROM构造了多种分析方法，并与标准算例进行了对比。虽然针对计算条件变化(如改变初始迎角[9]、CFD模型[10]、粘性[11]和结构参数[12]等)对二维翼型颤振边界的影响开展了一定的研究，但是在跨声速内，翼型厚度变化对颤振边界的影响则尚未给出研究规律。

本文基于CFD技术，采用系统辨识方法中的ARMA模型构造非定常气动力降阶模型[7]，耦合结构运动方程，建立一套频域/时域气动弹性ROM分析系统。然后采用该方法详细研究在翼型最大厚度所在位置保持不变时翼型厚度对颤振边界的影响。

1 气动弹性系统降阶模型

1.1 气动力降阶模型

基于ARMA模型的非定常气动力降阶模型是一种基于时间序列的降阶模型[3]，可以表述为如下方程：

其中，fa(k)为第k步得到的广义气动力；na、nb分别为系统输出和输入的延迟阶数；q表示输入的广义位移；Ai、Bi是模型中待辨识的参数。采用基于过滤高斯白噪声(FWGN)的模态激励信号进行模型辨识，该信号具有频带宽且幅值范围广的特点。

在构建气动力降阶模型时，采用MIMO (Multiple Input Multiple Output，MIMO) 系统方法，即一次输入n个位移信号同时激励n个模态的方式，通过最小二乘法拟合CFD求解出来的曲线，辨识出Ai、Bi矩阵。

定义状态向量xa(k)如下：

将离散差分方程转化为离散状态空间方程形式，则气动力降阶模型的状态空间方程和输出方程可以表示为：

其中，

1.2 气动弹性降阶模型

气动弹性系统的结构运动方程对应的状态空间方程可以表达为如下形式：

其中，

式中，q表示位移，M表示质量矩阵，K表示刚度矩阵，D表示阻尼矩阵。

将离散状态的气动力降阶模型转化为连续时间状态空间方程，并忽略静平衡气动力，即

式中，各字母含义与式(2)中一致，表示在连续时间状态下的气动系统矩阵系数。

耦合结构运动方程和气动力降阶模型方程，得到气动弹性系统的状态空间形式的控制方程：

其中，Q表示动压，状态空间矩阵如下：

可以对状态空间矩阵求解不同动压下的特征值，并采用模态跟踪技术[13]画出根轨迹，根据根轨迹分析系统的稳定性。

将上式转化为状态空间方程的形式通过时间推进求解[14]，如下：

其中，

采用叶正寅和张伟伟等人[15]构建的四阶隐式Adams线性多步法，通过预估-校正技术求解该隐式方程，提高气动弹性ROM的分析效率和鲁棒性[14]。

综上，建立基于系统辨识技术的气动弹性时域分析系统的耦合流程具体步骤如图1所示。

1.3 模型验证

为验证本文气动弹性ROM的精度，选取气动弹性分析的标准算例——Isogai案例A模型，该翼型由于在跨声速状态下呈现复杂的稳定性特性而被广泛应用于气动弹性程序验证。Isogai模型的翼型为NACA64A010翼型，是Isogai[16]提出计算颤振的二维翼型。计算中CFD求解器选取k-ω SST湍流模型，按照文献[17]，在所有马赫数下，基于弦长的雷诺数均给定为6×106。

在无粘条件下计算结果呈现为“S”形颤振边界，与参考文献结果一致。粘性颤振边界数值结果表明，跨声速凹坑最低点对应的马赫数大约为0.83，随后无因次颤振速度急剧增大；在Ma=0.87时无因次颤振速度开始减小，在Ma=0.9时达到第二个凹坑最低点，随后随着马赫数增加无因次颤振速度继续增大；当Ma≥0.93时，无因次颤振速度随马赫数的增加变化很小，基本保持在3.8左右。

图1基于系统辨识技术的气动弹性时域分析系统的耦合流程

上述结果是采用Intel(R) Xeon(R) CPU E5-2680 @2.50 GHz共19个进程并行计算，表1给出了直接采用CFD/CSD方法和基于ROM的频域分析方法的计算时间对比。可以得到，直接采用CFD/CSD方法计算耗时约7200 min，采用ROM预测颤振边界耗时246.4 min，计算效率提高1～2个数量级。

(a)无粘条件

(b)粘性条件

图2 采用CFD和ROM得到的Isogai翼型无因次颤振速度结果对比变化趋势

2 翼型厚度对颤振边界影响

翼型厚度影响着翼型的气动性能，是飞行器精细化设计时的重要参数，因此有必要研究翼型外形参数对颤振特性的影响。本小节研究了在最大厚度位置保持不变时，翼型厚度对颤振边界的影响。选取的翼型分别为NACA64A008、NACA64A009、NACA64A010、NACA64A011和NACA64A012，采用粘性条件，计算所用的气动参数和结构参数与1.3小节中的保持一致。

翼型厚度对颤振边界的影响如图3所示，图中展示了在最大厚度位置保持不变时，不同厚度翼型的无因次颤振速度边界和频率比边界。不同表面压力分布和摩阻系数分布对比图如图4～8所示。

数值结果表明：

(1)当马赫数Ma<0.80时，随着翼型厚度增加，颤振速度略有减小，但减小幅度在10 %以内。从图4可以看出，此时由于厚度的增加，导致激波强度逐渐增加，进而使得颤振速度逐渐减小。

(2)当马赫数Ma>0.91时，随着翼型厚度增加，颤振速度逐渐增加；且同一翼型随着马赫数的增加，颤振速度基本保持不变。从图8可以看出随着厚度的增加，分离区逐渐增加，激波强度基本保持不变，因此颤振速度逐渐增加。

(3)跨声速凹坑随着翼型厚度的增加而逐渐左移。从图5～7可以看出，在分离区作用下，跨声速凹坑左移，也就是第一次颤振速度增加。当进一步增加马赫数时，分离区域达到翼型后缘且随着马赫数增加而分离区减小，此时在激波的作用下，无因次颤振速度边界出现第二个台阶。

3 结语

本文采用ARMA模型构建了非定常气动力ROM，耦合结构运动方程，建立了频域/时域的气动弹性ROM。采用该方法研究Isogai翼型的跨声速颤振问题，并与CFD/CSD计算结果对比，结果表明该方法求解结果可靠。研究了最大厚度位置保持不变时，翼型厚度对颤振边界的影响，数值结果表明：当Ma<0.80时，随着翼型厚度增加，颤振速度略有减小；当马赫数Ma>0.91时，随着翼型的厚度增加，颤振速度逐渐增加。跨声速凹坑随着翼型厚度的增加而逐渐左移。

综上所述，采用ROM可以精确快速地预测颤振边界。当翼型最大厚度所在位置保持不变时，为了达到提高颤振速度的目标，可以尝试采取改变机翼翼型厚度的方式提高机翼对飞行环境的适应能力。