非定常气动力计算与颤振分析

苏 泽 杨红杰

（中航西安飞机工业集团股份有限公司西飞设计院，陕西 西安 710089）

在设计现代飞行器的过程中，必须要考虑气动弹性的影响。传统的气动弹性问题主要指发散、舵面反效代表的静气动弹性问题与颤振代表的动气动弹性问题。这类问题一般直接在频域中求解气动弹性方程，对亚音速飞行器来说，目前工程中普遍采用有限元方法直接进行模态分析，再利用偶极子格网法（DLM）求得频域离散形式的非定常气动力系数矩阵，从而实现气动力与结构模态间的直接耦合。

当飞行器包括控制系统时，必须要着重考虑控制系统与弹性机体结构之间的耦合作用，也就是气动伺服弹性力学（ASE）。传统飞控系统通常采用SISO 控制方式，工程中仍采用频域分析方法，基于经典控制理论进行稳定性分析。而现代飞行器控制回路之间的耦合作用明显，其飞控系统模型也具有MIMO 的特征，由此产生了以状态空间方法为基础的现代控制理论。此时，通过有理函数拟合方法将频域离散矩阵形式的非定常气动力转化为有理函数形式，使其能够进行拉氏反变换，从而建立状态方程，这是建立气动伺服弹性状态空间模型的关键[1]。

1 非定常气动力

通常情况下，气动弹性问题中的气动力主要是指飞行器弹性模态、舵面刚体偏转模态等广义位移带来的气动下洗产生的非定常气动力。

飞行器气动弹性一般运动方程如公式（1）所示。

式中：q为飞行器弹性模态坐标向量；δ为控制面刚体偏转坐标向量；ρ为大气密度；V为飞行速度；Mqq和Mqδ分别为飞行器弹性模态和舵面刚体偏转模态对应广义质量矩阵；Cqq为广义阻尼矩阵；Kqq为广义刚度矩阵；Qqq和Qqδ分别为弹性模态和舵面刚体偏转模态产生的广义非定常气动力系数矩阵。

对亚音速飞行器来说，工程中一般采用偶极子格网法来计算频域离散矩阵形式的非定常气动力。偶极子格网法是一种基于小扰动线化位势流方程的面元法，当计算非定常气动力时，需要合理地对气动面进行网格划分，将气动面沿展向和弦向分成若干个两侧边平行于来流方向的梯形面元，一般在面元网格的1/4 弦线上布置压力偶极子，在网格的3/4 处布置控制点，通过物面边界条件计算控制点的下洗速度[2]。

根据广义气动力的定义，其矩阵表达形式如公式（2）所示。

式中：R为广义气动力；NP为网格控制点处的模态矩阵，可以通过气动控制点与结构节点之间的插值矩阵，由结构模态矩阵转换得到；S为网格面积的加权矩阵，其对角项为各气动网格的面积；∆p为气动面元网格的压力分布。

根据非定常气动力理论，根据网格控制点满足的积分方程可以得到公式（3）。

式中：wq和wδ分别为弹性模态和舵面偏转模态带来的气动网格控制点的下洗角；D-1为气动力影响系数矩阵，也称AIC矩阵。

当计算飞行器弹性模态与舵面刚体偏转模态引起的非定常气动力系数时，气动网格控制点的下洗角与广义坐标之间满足公式（4）。

式中：b为参考半弦长；k为减缩频率，k=ωb/V；Nq和Nδ分别为网格控制点处的弹性模态和舵面偏转模态。

将公式（2）、公式（3）代入公式（4），得到广义非定常气动力的表达式，如公式（5）所示。

当给定气动模型的网格划分、飞行马赫数以及减缩频率时，可以确定飞行器弹性模态、舵面刚体偏转模态等带来的非定常气动力。

2 有理函数拟合

当采用时域方法分析气动伺服弹性问题和一些非线性问题时，要先对气动弹性广义运动方程进行拉氏变换，再将其转换为状态空间方程形式。而气弹运动方程中采用偶极子格网法求得的离散矩阵形式的频域气动力是减缩频率的复函数，无法直接进行拉氏反变换。此时，便需要引入有理函数拟合方法，在不损失计算精度的前提下，延拓频域气动力向拉氏域。

工程中最常用的有理函数近似方法包括ROGER 法、修正矩阵法（MMP）和最小状态法（MS）等，这些方法均以最小二乘法为基础。当采用这些方法拟合非定常气动力建立状态空间方程时，均需要引入气动力扩充项作为状态变量。不同方法引入的气动力扩充项数量也不同，ROGER 法产生的气动扩充项数量为结构模态数与气动滞后根的乘积，MMP法对应的气动力扩充项数量为广义气动力系数矩阵各列对应的气动滞后根数量之和，而MS 法对应的气动力扩充项数等于气动滞后根的数量[4]。研究表明，当气动力扩充项的数量相同时，MS 法的拟合精度最高，下面给出MS 法的主要拟合过程。

当给定飞行马赫数时，频域广义气动力系数矩阵可以为减缩频率的复函数，统一用Qq（k）表示。采用MS 法进行有理函数拟合的表达式如公式（6）所示。

将公式（6）展开为频域内实部和虚部的形式，如公式（7）所示。

当拟合精度较高时，选定的一系列减缩频率下的气动力拟合结果应与频域计算结果近似相等。为了获得更精确的有理函数拟合结果，需要假设在若干减缩频率点处的拟合结果与频域计算结果精确相等，即假设一定的约束条件：1）当约束减缩频率为0 时，气动力拟合结果与频域气动力相等。2）约束在减缩频率k1处气动力实部拟合结果与频域气动力实部相等。3）约束在减缩频率k2处气动力虚部拟合结果与频域气动力虚部相等。

由上述约束条件及其他约束条件可以将拟合公式中A0、A1和A2分别转化为减缩频率k、D和E矩阵的表达式，以进一步给出其他非约束减缩频率点处的气动力系数拟合公式。

为了进一步求得拟合公式的解，MS 法需要先给定R的矩阵元素，再由最小二乘法确定矩阵D和E。首先，给定矩阵E，按行拟合出矩阵D。其次，由现有的矩阵R和D，按列拟合求出矩阵E。最后，计算拟合的精度，如果拟合的精度不满足要求，就重复前面的拟合过程，反复迭代计算D-E-D，直到得到满意的拟合结果。一般情况下，迭代10次即可收敛。

当求出各系数矩阵后，令s=ik（s为拉普拉斯算子），将减缩频率转化为拉氏变量，气动力拟合如公式（8）所示。

由此，利用有理函数拟合方法就可以将广义非定常气动力从坐标轴拓展到整个拉氏域。

3 算例

该文以一个典型翼面为例对上述非定常气动力计算方法和有理函数拟合方法进行验证。

3.1 计算模型

计算模型由机翼翼面与机翼根部连接组成，翼面采用板单元建模，根部连接为梁单元，并在梁单元根部固支约束。利用MSC.NASTRAN 软件进行模态分析，得出前五阶结构弹性模态见表1。

表1 模态分析结果

采用ZAERO 软件建立机翼翼面的气动网格，为了简化步骤，沿弦向和展向均等分为5 段。为了验证有理函数拟合方法的精确性，分别采用频域法和状态空间方法对上述模型进行开环颤振分析，颤振计算方法选取ZAERO 软件中的g法，固定高度为海平面，马赫数0.05，计算选取前五阶弹性模态，不考虑阻尼，取减缩频率为0.0～1.0。

当状态空间法计算时，采用MS 法进行广义气动力系数矩阵的有理函数拟合，当约束减缩频率为0 时，气动力拟合结果与频域气动力相等，取拟合迭代次数为100 次，并根据经验公式法选取10 个气动滞后根，如公式（9）所示。

式中：Ri为第i个气动滞后根取值；kmax为计算选取的最大减缩频率；Nlag为气动滞后根数量。

3.2 颤振结果对比

采用ZAERO 软件中的g 法进行颤振分析，得到频域法V-G、V-F 曲线，如图1所示。2 种方法颤振频率与颤振速度对比见表2。

由图1 可知，计算得到一支颤振，其类型为典型弯扭耦合颤振，参与模态为机翼一阶弯曲和一阶扭转模态，颤振频率为8.112 Hz，颤振速度较高。在采用MS 法进行有理函数拟合后，颤振频率及颤振速度与频域计算结果偏差较小，说明气动力拟合过程并没有降低计算精度。

3.3 非定常气动力拟合对比

频域法计算不采用有理函数拟合，其得到的非定常气动力系数为精确解，可为n×m维离散矩阵，每项矩阵元素均包括实部和虚部2 个部分（n为计算采用的弹性模态数量；m为弹性模态数量与舵面刚体偏转模态数量的和）。该算例中不包括操纵面，机体弹性模态数量为5，因此频域法计算得到的非定常气动力系数矩阵为5×5 的复数矩阵。

第3.2 节得出颤振的主要参与模态为翼面的第一阶模态、二阶模态，不同减缩频率的2 种方法得到的广义气动力系数矩阵元素Q11与Q22的对比如图2所示，下标表示元素在矩阵中的位置。

图2 矩阵元素随减缩频率变化趋势对比

对比气动力系数矩阵元素可以得出，在给定的减缩频率下，采用MS 法进行非定常气动力有理函数拟合，当选择合适数量的气动滞后根时，拟合得到的广义气动力系数矩阵元素与频域直接计算值已基本一致，气动力随减缩频率的变化规律也与频域法一致，拟合过程整体偏差较小，拟合效果较好，结论与颤振对比结论一致。

4 结语

该文论述了亚音速飞行器非定常气动力建模方法，给出了基于MS 法进行广义气动力有理函数拟合的计算方法，并根据典型翼面算例分析了频域气动力、拟合后的气动力对颤振结果的影响。分析结果表明，当采用MS 法进行非定常气动力拟合时，选择合适数量的气动滞后根后能够得到较准确的气动力计算结果，其拟合精度能够满足颤振分析需求。

实际上，针对线性系统的颤振分析、频域ASE 分析、机动载荷分析以及离散阵风分析等传统气动弹性问题，工程中一般直接采用频域计算得到的非定常气动力；而针对时域ASE 分析及其他非线性动力学响应分析，例如非线性颤振、投放载荷分析等，难以直接使用频域气动力，此时，必须要进行有理函数拟合，通过状态空间方法得到可靠结果。