非稳定动态过程非定常气动力建模

2020-09-10陈森林高正红朱新奇庞超杜一鸣陈树生

航空学报 2020年8期

陈森林，高正红，朱新奇，庞超，杜一鸣，陈树生

西北工业大学航空学院，西安 710072

对先进战斗机而言，在大迎角范围内能够进行大幅度可控的快速机动是一项必备的能力。当战斗机进行过失速机动时，绕飞机的流场会出现分离涡等复杂流动现象，此时作用在飞机上的气动力(矩)将呈现明显的非线性和非定常特性，即气动力(矩)不仅取决于当前时刻的运动状态，还与运动的时间历程有关[1]。对于这种情况，传统的气动导数模型已经很难适用[2]，需要探索更加准确高效的非线性非定常气动力建模方法。

相关建模方法已经有不少研究[3-5]，常见的包括阶跃响应模型[1]、Volterra级数模型[6-7]、状态空间模型[8-9]、微分方程模型[10-12]等。随着机器学习的蓬勃发展，近年来出现了将机器学习应用到气动力建模的尝试，如神经网络[13]、模糊逻辑[14]、支持向量机[15-16]等，这类方法的优点在于将系统视为黑箱，利用机器学习对复杂系统的描述能力预测输出。其中，支持向量机(Support Vector Machine, SVM)是Vapnik[17]在1995年基于统计学习理论提出的一种针对小样本情况的机器学习方法，其最初应用于模式识别，后来逐渐扩展到系统辨识、智能控制等领域[18-19]。最小二乘支持向量机(Least Squares Support Vector Machine, LS-SVM)[20]是标准SVM的一种改进方法，其将二次规划问题转化为线性方程组求解，降低了计算复杂度。

在现有的非定常气动力建模过程中，大多是以模型在风洞中的强迫振动试验所得到的数据为初始样本进行，其流程如图1所示。在某一减缩频率下，当试验对象经过多个周期的强迫振动之后，试验对象和其周围的流动以及产生的气动力均呈现出稳定的周期性变化特征，此时的运动是一个稳定的动态过程，此过程的气动力数据即是该减缩频率下的稳定滞环数据，采集该数据作为样本。接下来，在另一减缩频率下继续进行试验。当所有试验完成之后，获得多个减缩频率下的稳定滞环数据。最后，基于这些数据进行建模。建模完成之后，模型被用来预测其他减缩频率下的气动力响应。但是，这些方法针对的对象始终是稳定的滞环，对于气动力是如何达到稳定滞环的，也就是气动力的初始发展过程，即图1中气动力稳定之前的阶段没有考虑。飞机的快速机动过程是一个动态的变化过程，可以认为在这个过程中飞机不会持续运动到进入稳定滞环，而是处于初始的非稳定阶段。因此，初始非稳定阶段的气动力预测相比稳定段，更具有实际意义。这里以某机翼为例，由CFD计算的在不同减缩频率k下做俯仰运动时的升力系数CL滞环，如图2所示，初始迎角α为30°。图中，实线表示稳定滞环，虚线表示进入稳定滞环的初始非稳定过程。可以发现，气动力明显存在一个进入稳定滞环的过程。所以，仅仅预测稳定滞环是不够的，更重要的是能够预测初始非稳定过程的气动力。

图1 基于风洞试验数据的气动力建模方法示意图

图2 某机翼做俯仰运动时的升力系数滞环(Ma=0.4,Re=3.0×106)

传统方法基于多个减缩频率的稳定滞环数据进行建模，存在的另一个问题是，系统每次仅受到单个频率的激励，只能激发出系统的局部线性特性，而大迎角范围内的气动力是非线性的。对于战斗机的过失速机动飞行而言，真实的飞行过程以纵向俯仰为例，战斗机在拉升的过程中快速进入过失速机动，该过程几乎不可能是标准的振动，更不可能用单一频率的稳定振动进行模拟。为了能够正确地模拟战斗机动态进入过失速机动过程以及产生的非定常气动力(矩)，不仅需要能够包含非稳定动态特性影响，同时需要反映不同频率的非线性耦合影响。

本文根据Morelli等[21-22]考虑多操纵面非线性耦合影响而提出的正交相位优化多正弦激励方法，并合理选择输入形式，设计了用于构建大迎角俯仰动态过程的非线性激励模型。以某机翼为例，通过CFD技术获取了激励模型下的气动响应，并以此为构建大迎角非定常气动力模型的样本数据，利用LS-SVM方法建立了仅通过一次激励，就能预测在幅值和频率范围内任意运动的非线性非定常气动力模型。算例结果表明，该模型不仅能够正确地预测不同频率下的稳定滞环，同时能预测初始的非稳定过程。

1 激励模型

1.1 非稳定动态过程和非线性系统辨识

首先借助振动理论讨论一般线性系统在外部激励下的响应，分析非定常气动力响应的变化过程。以一个典型的二阶振动器系统为例，其运动方程为

(1)

式中：y为位移；ζ为阻尼比；p为自然频率；Y0p2cos(ωt)为外部输入。

基于振动理论，位移响应y包含两部分，分别为自由振动y1和强迫振动y2，即

y=y1+y2

(2)

y1=e-ζpt(Asin(qt)+Bcos(qt))

(3)

y2=Ycos(ωt-φ)

(4)

对于一个稳定系统，自由振动y1会随时间的推移而逐渐衰减到零，最后仅剩下强迫振动y2。以某个振荡器[23]为例，不考虑非线性环节时，在外部输入为正弦的情况下，系统的自由振动响应和强迫振动响应如图3所示。可见，系统响应在开始的时候，存在一个由自由振动响应和强迫振动响应共同构成的非稳定段，这个过程是一个非稳定动态过程。随着时间的推移，自由振动衰减到零，只剩下强迫振动，此时的响应达到稳定，这个过程是一个稳定动态过程。传统的非定常气动力建模方法，仅采用了稳定段的非定常气动力数据，而忽略了初始的非稳定过程。换言之，只有强迫振动数据被采集和使用，而自由振动数据被遗漏了，这将导致对系统响应变化过程的信息缺失。因此，基于这些数据的模型只能预测最终的稳定滞环，而不能预测进入滞环的初始非稳定过程。

图3 某振荡器[23]在正弦激励下的强迫振动响应和自由振动响应

其次，大迎角区的气动力会呈现明显的非线性特征。对于非线性系统而言，单个频率激励只能激发出系统在该频率下的局部特性，无法激发出系统中不同频率耦合作用的影响。因此，如果仅分别使用不同频率的稳定滞环为样本数据进行大迎角非定常气动力建模，所获得的模型只能够反映不同频率下稳定的非定常气动力的迟滞环特性，不能反映系统中频率耦合的非线性影响。

因此，为了建立能够正确计算过失速机动过程的非定常气动力模型，首先需要构建能够激发非线性动态过程中频率耦合作用的激励模型，其既包含了非稳定动态过程的影响，又适于非线性系统辨识。

1.2 激励的选择和设计

对于系统辨识和机器学习问题，激励的选择和设计至关重要，直接影响了模型的预测精度。从系统辨识的角度，选择训练数据的关键是充分激发出系统特性。对于线性系统辨识，有多种输入，例如扫频[24]、3211[25]等可以采用。对于非线性系统辨识，缺乏通用的输入信号。

根据Prazenica和Kurdila[26]的研究，在频域表示线性系统特性的频率响应函数是一阶的，可以表示为H1(jω)，其中ω是频率，j是虚部。以扫频激励为例，其触发到了不同频率下的系统特性，因此适用于辨识一阶频率响应函数H1(jω)，即适用于辨识线性系统。表示非线性系统特性的频率响应函数是高阶的，包含了H1(jω)、H2(jω1,jω2)、H3(jω1,jω2,jω3)等。换言之，非线性系统的响应不仅包含了单一频率的影响，还包含了不同频率的耦合影响，此时扫频激励不再适用。因此，激励的选择和设计必须能够反映不同频率的耦合作用。

这里采用Morelli等[21-22]设计的用于飞机多操纵面耦合激励的输入——正交相位优化多正弦，作为训练输入。为了对飞机多个操纵面同时进行激励，以实现飞机多通道气动导数的同时辨识，Morelli等需要使得不同操纵面的激励在时域和频域上都互相正交，以实现互不干扰。这种输入通过在不同通道将多个互相正交的正弦信号叠加，然后利用优化方法调节相位以避免不同分量相加使得幅值过大，导致系统受到激励后偏离参考点过远，最终使得不同通道之间的激励互相正交，同一通道中不同分量之间也互相正交，并覆盖感兴趣的频段。借鉴这种思路，这里将多通道合并到一个通道，同时仍然保持正交性，构造了一种能体现不同频率耦合作用的输入，用于非线性系统辨识。本文作者[27]成功将这种输入应用到Volterra级数高阶核的辨识，实现了翼型在跨声速下小迎角范围内的非线性非定常气动力建模。

正交相位优化多正弦输入的表达式为

(5)

式中：u为输入；Ai为幅值；T为输入时长；φi为相位；M为选择的正弦总数。

根据Schroeder[28]的研究，激励信号要尽可能充分的激励系统，需要激励信号的能量尽可能大，衡量信号能量的是均方根(Root Mean Square, RMS)，表示为

(6)

式中：N为输入的离散点数。同时，还要使得系统响应不偏离参考点过远，保证这一点需要激励的幅值差值尽可能小，即max(u)-min(u)尽可能小。综合考虑这两个因素，Schroeder提出，可以利用相对峰值因子(Relative Peak Factor, RPF)最小来确定输入相位，其定义为

(7)

该输入设计的具体过程如下：

步骤1根据研究对象的实际情况确定输入的时长T和频带[fmin,fmax]，一般要求fmin≥2/T，由此可以得到频率精度Δf=1/T和M=fix[(fmax-fmin)/Δf]+1，其中fix(·)表示向下取整。

步骤3随机给定各个分量的相位初值φi∈(-π,π]，然后运用优化方法调节相位，使得各分量相加后的输入的相对峰值因子最小，最后得到需要的相位分量φi。

步骤4根据前3步确定的参数构成的输入作小幅相位偏移，使得输入从零开始，到零终止，便于对系统施加激励，最终得到需要的输入。

图4 正交相位优化多正弦输入

1.3 训练所用输入形式的选择

在训练之前，必须根据实际问题的物理背景选择合理的输入形式。因为非定常气动力不仅依赖于当前的运动状态，还依赖于之前的运动过程，所以输入形式的选择必须反映这一点。在已有的利用机器学习进行气动力建模的研究中，减缩频率或者人为定义的等效减缩频率常常作为输入参数之一。例如，Ignatyev[29]、Wang[15]和史志伟[30]等，所用训练数据均来自于不同减缩频率下的强迫振动试验。然而，如前所述，真实的飞行过程几乎不可能是标准的振动，因而不存在固定频率。建立一个不依赖于频率且具有通用性的气动力模型是本文研究的问题之一。

这里选择的输入形式为

x=[α(τ),α(τ-1),…,α(τ-m)]T

(8)

式中：α为迎角；m为当前时刻的气动力对之前运动过程的依赖长度。此处没有引入减缩频率，因为在激励信号的设计中已经通过频带的选择考虑了频率的影响。

2 最小二乘支持向量机

2.1 基本理论

给定训练样本集(xi,yi)，i=1,2,…,n，xi∈Rd，yi∈R，SVM回归的基本思想是将样本所在的输入空间Rd通过非线性映射φ映射到特征空间φ(x)，然后在该空间构造回归函数y=f(x)wTφ(x)+b，其中w为加权系数，b为偏差。参数(w,b)通过结构风险最小化原则确定，即求解如下优化问题：

(9)

用Lagrangian方法求解式(9)的等式约束优化问题，得

(10)

式中：λi为Lagrangian乘子。最优解为

(11)

对应条件为

(12)

(13)

cξi=λi

(14)

yi=wTφ(xi)+b+ξi

(15)

将式(12)和式(14)代入式(15)，消去w和ξi，得

(16)

式中：〈· 〉为内积。定义核函数K(xi,xl)=〈φ(xi),φ(xl)〉，将式(13)和式(16)合并，得方程组的矩阵形式为

(17)

式中：Kil=K(xi,xl)；λ=[λ1,λ2,…,λn]T；y=[y1,y2,…,yn]T；1n为n个1的列向量。

为了求解式(17)以得到参数(λ,b)，常用方法是利用核函数矩阵K的对称正定性质对方程组做变换，然后利用Cholesky分解求解，具体方法可见文献[31]。已知(λ,b)，即得到了所需的LS-SVM模型：

(18)

2.2 模型参数的确定

SVM通过引入核函数避免了在高维特征空间直接计算内积，极大简化了计算。根据泛函的有关理论，任意满足Mercer条件的对称函数都能作为核函数，最常采用的是高斯核函数：

(19)

核函数选定后，模型参数的确定主要是核参数σ和正则化参数c的确定。某组参数(σ,c)所对应的泛化误差常采用留一交叉验证法(Leave-One-Out Cross-Validation)计算，即给定一组(σ,c)之后，每次将训练样本中的一个样本作为测试集，剩下的作为训练集，计算预测误差，如此重复，直到每个样本都作为过测试集，最后将所有预测误差的平均值作为该组参数下的泛化误差。如何找到使得泛化误差最小的参数(σ,c)，可采用MATLAB优化工具箱中的fmincon函数，其利用共轭梯度方法进行搜索[31]。

2.3 训练方法

输入信号选择以后，接下来利用CFD计算气动力。然后根据选择的输入形式构造输入向量xi，输出分别为升力系数、阻力系数和俯仰力矩系数。最后对LS-SVM进行训练。具体步骤如下：

步骤1按照1.2节构造输入，然后由CFD计算气动力响应。

步骤2根据1.3节构造输入向量xi，输出yi为气动力响应，i=1,2,…,n。

步骤3给定模型参数的初值，包括式(9)中的c和式(19)中的σ。

步骤4利用MATLAB的fmincon函数，以留一交叉验证法计算的平均预测误差为优化目标，寻找最优的模型参数(σ,c)。

步骤5回到步骤3和步骤4，可重复多次，避免(σ,c)的局部最优，最终选择使得平均预测误差最小的模型参数(σ,c)。

步骤6求解式(17)，得到参数(λ,b)，即完成了模型的训练。在已知测试数据后，由式(18)即可计算模型的预测结果。

3 训练

以某机翼在亚声速下大迎角范围内做俯仰运动为例，验证LS-SVM模型预测非线性非定常气动力的能力和所采用设计输入的有效性。训练所用数据由CFD计算得到，与所用模型无关。计算条件为Ma=0.4，Re=3.0×106，采用刚性动网格。机翼外形参数如图5所示，剖面翼型为NACA-64A204，网格图如图6所示。

图5 机翼外形图

图6 机翼网格生成图

本文所用CFD求解器是NASA的Langley研究中心开发的开源计算工具CFL3D[32]，其采用有限体积法和Spalart-Allmaras湍流模型求解雷诺平均Navier-Stokes (RANS)方程。CFL3D已经在航空领域通过多个算例得到了广泛验证[33-35]。由CFL3D计算的机翼在当前状态下的定常气动力如图7所示。

图7 机翼在Ma=0.4、Re=3.0×106下的定常气动力系数

以设计的输入由CFD进行非定常计算，得到气动力数据，然后进行训练。训练完成后，以训练数据作为测试样本检验训练效果，图8和图9分别为以不同基准迎角做俯仰运动下，升力系数、阻力系数和俯仰力矩系数的LS-SVM预测结果和CFD计算结果。Case 1的基准迎角为0°，幅值为80°；Case 2的基准迎角为30°，幅值为30°。关于式(8)中时间依赖长度m的选择，需要尝试不同的时间长度，然后根据训练误差进行选择，这里Case 1取m=10，Case 2取m=30。由图8和图9可见预测结果和CFD结果非常吻合，验证了LS-SVM对训练数据的学习能力。

图8 Case 1：初始迎角α0=0°时LS-SVM预测结果与训练数据对比(前10 s)

图9 Case 2：初始迎角α0=30°时LS-SVM预测结果与训练数据对比(前10 s)

4 非定常气动力验证

模型训练完成后，需要验证其在其他输入下的预测准确性，也即验证其推广能力。这里分别对机翼在不同基准状态下俯仰做正弦和扫频运动时的气动力进行计算验证。

4.1 Case 1

Case 1的基准迎角为0°，幅值为80°，机翼做俯仰运动时迎角变化为

(20)

前3种正弦运动对应的减缩频率k分别为0.01、0.03、0.06，第4种为扫频运动，对应的减缩频率从0线性增加到0.08。

由LS-SVM计算的测试运动的气动力响应和CFD验证结果对比如图10所示。3种正弦运动的气动力滞环包括初始非稳定动态过程，都和CFD结果一致。图中虚线加粗部分是进入滞环的非稳定动态过程，可以发现，当减缩频率为0.01和 0.03时，气动力的非稳定段与稳定段几乎重合，差异不明显；当减缩频率增加到0.06时，气动力的非稳定段与稳定段出现了差异。通过对比可以看到，不仅稳定滞环部分得到了准确预测，进入稳定滞环的初始非稳定过程也得到了准确预测。

图10 初始迎角α0=0°时LS-SVM预测结果与CFD结果

除了3种正弦运动，扫频运动的预测结果与CFD结果也吻合很好，表明通过采用文中的训练数据设计方法和输入形式，避免了以减缩频率为参数，不仅能预测正弦运动，理论上能预测范围内的任意运动，使得模型更具有普适性。

4.2 Case 2

Case 2的基准迎角为30°，幅值为30°，机翼做俯仰运动时迎角变化为

(21)

前3种正弦运动的减缩频率分别为0.01、0.03 和0.06，第4种为扫频运动，对应的减缩频率从0线性增加到0.08。

由LS-SVM计算的测试运动的气动力响应和CFD验证结果对比如图11所示。3种正弦运动中的气动力滞环随着减缩频率的增大，非稳定动态过程表现逐渐明显。包括非稳定动态过程在内的整个滞环，都实现了准确预测，即使升力系数、阻力系数和俯仰力矩系数随减缩频率的变化呈现不同的特征，预测结果与CFD结果都非常吻合。

图11 初始迎角α0=30°时LS-SVM预测结果与CFD结果

与Case 1类似，扫频运动的预测结果与CFD结果也吻合很好，再次表明通过采用文中的训练数据设计方法和输入的形式选择，不仅能预测正弦运动，理论上能预测幅值和频率范围内的任意运动，使得模型更具有普适性。

Case 2相比Case 1，基准迎角发生了变化，气动力的初始非稳定过程表现得更为明显，也就是说，在Case 2状态下，气动力的初始段和稳定滞环存在明显差异。根据前面的分析，飞机的真实机动过程始终处于初始段，不会到达稳定滞环。因此可以认为，在当前状态下，仅依靠稳定滞环来预测飞机机动过程中的非定常气动力将存在很大误差。

4.3 以强迫振动为训练输入时的预测

现在将1.1节讨论的强迫振动作为训练输入，其他不变，选用式(21)中的第2个正弦运动为算例，验证气动力的预测效果。由CFD计算的升力系数的前3个周期响应如图12所示。可以看到，气动力很快趋于周期性变化。因此，可以认为后两个周期的气动力响应是稳定的。分别将后两个周期的稳定振动数据和包含初始非稳定段在内的全部3个周期的振动数据作为训练数据，然后预测同样正弦运动下的升力系数响应。预测结果的时间响应曲线(前0.6 s)和滞环曲线如图13所示。标注“稳定输入”的曲线为采用后两个周期的稳定数据为训练数据所得到的预测结果，标注“全部输入”的曲线为采用全部3个周期的数据为训练数据所得到的预测结果。

图12 初始迎角α0=30°时由CFD计算的升力系数训练输入(减缩频率k=0.03)

图13 初始迎角α0=30°时不同训练输入下LS-SVM预测结果与CFD结果(减缩频率k=0.03)

通过和CFD结果对比可以发现，采用稳定振动数据进行训练，可以准确预测稳定滞环，但进入稳定滞环的初始非稳定过程无法实现预测。采用包含初始变化过程的振动数据进行训练之后，可以实现气动力整个发展过程的准确预测。这个例子再次表明，采用强迫振动输入只能预测到稳定滞环，只有当包含初始的自由振动之后，进入滞环的非稳定过程才能得到预测。因此，基于强迫振动试验数据的建模方法难以用于飞机快速机动过程中非定常气动力的预测。