飞行安全是大型商用运输机设计研究中备受关注的核心问题。按照国际民航组织(ICAO)和商业航空安全委员会(CAST)分类标准^[1]，飞行事故按类型分为：可控飞行撞地(CFIT)、空中失火(F-NI)、燃油泄漏(FUEL)、飞行失控(Loss of Control-in Flight, LOC-I)、空域冲突(MAC)、冲出跑道(RE)、动力系统故障(SCF-PP)等。按事故类型统计商用航空运输中人员伤亡飞行事故显示^[2]，飞行失控在2008—2017的10年间共造成14起事故、1 131人死亡，在所有飞行事故类型中，无论事故数量还是死亡人数均占比最高。对飞行失控事故的飞行记录数据分析表明^[3]，因恶劣天气、系统故障或机组人为因素等使飞机进入超出正常飞行包线的极限飞行状态，是造成飞行失控的重要原因之一。近年发生的波音737MAX事故进一步说明，在一些特定的民用飞机设计约束之下，飞机迎角的快速增加，并进而进入失速或过失速状态，是民用飞机在节能、减排等多种约束目标优化的情况下可能会出现的问题，而目前波音737MAX的处理方法存在很大的风险，因此需要从飞机本体特性和新飞行控制律的角度加以研究。事实上，按线性、定常条件设计的飞机遭遇非线性、非定常的动态气动力环境, 一直存在潜在的风险，如果对其动态气动特性缺乏足够的认识，没有建立精确的动态气动力模型以支撑驾驶员进行应对此种意外情况的应急培训，其结果往往是灾难性的。通过开展民机极限飞行状态的动态气动力特性研究, 来改善飞机飞行失控预防、极限状态改出、飞行模拟训练和飞行事故分析等, 以减少或杜绝类似事件再次发生, 无疑是今后减少和避免国内外诸多空难事故需要探索的一个重要技术途径。

飞机不慎进入超出正常飞机包线的极限飞行状态时，有两个显著特征^[4]：一是迎角、侧滑角变化范围广；二是飞机处于运动快速变化状态。这种复杂的非常规运动环境与人们熟知的定常、直线飞行有着本质区别，这不仅体现在飞机的气动力特性本身——气动力具有完全的动态特征，即非线性、非定常性质；而且体现在气动力风洞试验设计与动态气动特性研究的关系方面，即这种情况下的动态气动特性研究究竟需要怎样的风洞试验数据，以及风洞试验数据如何用于构建气动力模型的问题^[5]。近30年来，关于飞机的大迎角与动态气动特性研究，主要围绕战斗机大迎角过失速机动的气动力试验与建模开展，发展起来的动导数、大振幅振荡、旋转天平等动态气动力风洞试验方法^[6-7]，已经在战斗机过失速机动与尾旋研究中得到广泛应用。相应的，发展了多种建模方法，包括以气动导数模型^[8]、非线性阶跃响应模型^[9]、状态空间模型^[10]、微分方程模型^[11]等为代表的数学建模方法以及以神经网络模型^[12]、模糊逻辑模型^[13]和支持向量机模型^[14]等为代表的人工智能建模方法，这些方法应用在战斗机布局飞机气动特性研究与建模中，使得战斗机大迎角机动过程中的动态气动力问题得到相当充分的认识和不同程度的解决^[15-17]。但是，动态气动力特性与飞机布局形式密切关联^[18]，而目前对于运输类布局飞机极限飞行状态气动力试验与建模研究仍然有限，对于民机失速，尤其是过失速动态气动特性尚缺乏深入的风洞试验研究；目前飞行模拟器中对于超出正常迎角/侧滑角包线范围的数据，一般是在正常包线数据基础上外推^[19]，存在定性的误差，也难以支撑应对极限飞行状态的飞行模拟训练。

本文针对典型民机布局飞机模型，主要关注极限飞行状态中动态气动力的非线性与非定常特征，从风洞动态试验方法中蕴含的假设和应用条件出发，结合极限飞行状态特点和动导数、大振幅试验结果，分析大迎角区域动态气动力参数影响规律和非线性模型结构，在此基础上探讨极限飞行状态动态气动力的风洞试验设计和非定常气动力建模问题。

1 飞机模型

考虑到布局典型性，选择NASA通用运输机标模(Common Research Model, CRM)作为研究模型。CRM是NASA发布的代表典型双发、远程、双通道宽体商用运输机布局的标模，飞机三维数模及数据面向国际合作公开发布^[20]。该标模对机翼进行全新设计，采用现代先进的超临界翼型，飞机机翼展弦比为9.0，根稍比为0.275，机翼1/4弦线后掠角为35°；而飞机的机身、平尾和垂尾等部位的关键尺寸、布局参数与波音777-200保持一致^[21]，设计巡航马赫数Ma=0.85(设计点升力系数C_L=0.5)，飞机气动布局如图 1所示。值得一提的是，在NASA发布的CRM原始标模中，飞机不带操纵面，本研究中为了后续开展极限飞行状态操控特性研究需要，参考波音777-200操纵面设计了升降舵、副翼和方向舵。

采用2.45%缩比模型开展动态气动力风洞试验，模型主要参数如表 1所示^[22]。

2 试验方法 2.1 动态气动力数学模型及其简化形式

飞机受到的气动力可以严格地表示为飞行状态变量的泛函：

$ \begin{array}{*{20}{c}} {{C_i}(t) = {C_i}(M(t),h(t),\delta (t),\alpha (\tau ),\beta (\tau ),}\\ {p(\tau ),q(\tau ),r(\tau ))\;\: - \infty \le \tau \le t} \end{array} $	（1）

式中：i=L, D, Y, M_x, M_y, M_z，分别表示升力、阻力、侧向力和三轴气动力矩；C_i(t)表示在给定t时刻的气动力和力矩，其不仅与t时刻的飞行状态参数有关，而且与飞行状态参数变化的历史有关。这个表达式能够反映空气流动产生的气动力的非定常、非线性实质，但显然并没有简单的途径可以得到这种泛函的数值^[23]。为此，目前的工程实践中通常引入3个假设对这一泛函做如下准定常化、线性化处理^{[7, 23]}：

假设1   假设较远历史状态对t时刻的气动力影响可以忽略，以飞行状态参数α(τ)为例，在t时刻附近展开为泰勒级数：

$ \begin{array}{*{20}{l}} {\alpha (\tau ) = \alpha (t) + (\tau - t)\dot \alpha (t) + }\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \frac{1}{2}{{(\tau - t)}^2}\ddot \alpha (t) + \cdots } \end{array} $

（2）

则可以用t时刻的α及其各阶导数值代替式(1)中的α(τ)，其他飞行状态参数作类似处理:

$ \begin{array}{*{20}{l}} {{C_i}(t) = {C_i}(M,h,\delta ,\alpha ,\dot \alpha , \cdots ,\beta ,\dot \beta , \cdots ,}\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} p,\dot p, \cdots ,q,\dot q, \cdots ,r,\dot r, \cdots )} \end{array} $

（3）

引入该假设后，气动力变成关于飞行状态参数的准定常、非线性模型，大部分实际问题中，只保留到飞行姿态角速度，就可以保证足够的精度，即

$ {C_i}(t) = {C_i}(M,h,\delta ,\alpha ,\beta ,\dot \alpha ,\dot \beta ,p,q,r) $	（4）

不失一般性，式(4)可以分解为两部分，即静态气动力C_i1和动态气动力C_i2，C_i1是关于飞机平动运动参数的非线性函数(与飞机姿态角相关)，C_i2是关于飞机转动参数的非线性函数(与飞机姿态角速度相关)：

$ \begin{array}{l} {C_i}(t) = {C_{i1}}(M,h,\delta ,\alpha ,\beta ) + \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {C_{i2}}(M,h,\delta ,\alpha ,\beta ;\dot \alpha ,\dot \beta ,p,q,r) \end{array} $	（5）

假设2   假设动态气动力中，各运动参数对气动力影响可以相互解耦，即

$ \begin{array}{l} {C_{i2}}(M,h,\delta ,\alpha ,\beta ;\dot \alpha ,\dot \beta ,p,q,r) = \\ \begin{array}{*{20}{l}} {{C_{i2}}(M,h,\delta ,\alpha ,\beta ;\dot \alpha ) + {C_{i2}}(M,h,\delta ,\alpha ,\beta ;\dot \beta ) + }\\ {{C_{i2}}(M,h,\delta ,\alpha ,\beta ;p) + \cdots } \end{array} \end{array} $	（6）

假设3   假设动态气动力中，气动力随飞机姿态角速度是线性变化的，即

$ {C_{i2}}(M,h,\delta ,\alpha ,\beta ;p) = {C_{i\bar p}}(M,h,\delta ,\alpha ,\beta )\bar p $	（7）

式中：${C_{i\bar p}} = \partial {C_{i2}}(M, h, \delta , \alpha , \beta ;p)/\partial \bar p, \bar p = pb/2V$为无因次角速率。

由此，可得

$ \begin{array}{l} {C_{i2}}(M,h,\delta ,\alpha ,\beta ;\dot \alpha ,\dot \beta ,p,q,r) = \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {C_{i\dot \alpha }}\frac{{\dot \alpha c}}{{2V}} + {C_{i\dot \beta }}\frac{{\dot \beta b}}{{2V}} + {C_{ip}}\frac{{pb}}{{2V}} + {C_{iq}}\frac{{qc}}{{2V}} + {C_{ir}}\frac{{rb}}{{2V}} \end{array} $

（8）

综上，准定常条件下气动力线性化模型表达式为

$ \begin{array}{*{20}{l}} {{C_i}(t) = {C_{i1}}(M,h,\delta ,\alpha ,\beta ) + {C_{i\dot \alpha }}\frac{{\dot \alpha c}}{{2V}} + {C_{i\dot \beta }}\frac{{\dot \beta b}}{{2V}} + }\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {C_{ip}}\frac{{pb}}{{2V}} + {C_{iq}}\frac{{qc}}{{2V}} + {C_{ir}}\frac{{rb}}{{2V}}} \end{array} $

（9）

基于上述假设的数学模型以及发展的与之相应的风洞动态试验方法，一般都严格限制在小到中等迎角时气动力非线性、非定常现象很弱的范围，对于极限飞行状态，上述假设是否仍然适用及其所带来的影响，以及如何进行试验设计和建模需要结合特定飞机布局特点进行分析。

2.2 试验内容与方法

结合2.1节分析，尽管在大迎角条件下，上述准定常化、线性化方法，难以严格描述复杂非线性特征，但是实际上任何的非线性气动力现象都是从线性气动力出现异常发展而来的^[23]。基于这个认识，针对极限飞行状态动态气动力试验与建模，本文基本思路是：利用已经发展成熟的动导数试验方法，开展包含大迎角在内的不同频率、不同振幅的动导数试验，分析不同迎角区域参数影响规律和动态气动力的非线性特征，建立非线性动态气动力模型；利用大振幅试验结果，对比线性和非线性气动力模型，对动态气动力的非线性特征和非线性建模结果进行分析和验证，并进一步研究气动力的非定常特征，进行非定常气动力建模。

动导数试验通过强迫振荡方式进行，分为俯仰强迫振荡、偏航强迫振荡、滚转强迫振荡试验。试验在中国空气动力研究与发展中心(CARDC)FL-14风洞开展(图 2)。试验时，模型绕相应体轴在不同迎角作不同振幅、频率的正弦振荡，按照飞机的短周期或荷兰滚运动模态特征，或者飞机飞行运动中涉及的飞行状态参数典型范围，确定参数模拟范围。试验风速为30 m/s，试验迎角范围为-10°~60°，振荡频率范围为0.5~1.5 Hz，振幅范围为3°~10°。

大振幅振荡试验在同一套试验装置上开展，在不同迎角区域作振幅分别为10°、15°、20°的运动，频率范围在0.25~1.25 Hz之间。

3 动态气动力非线性特性 3.1 非线性特征分析

以通过俯仰振荡获取俯仰动导数为例进行分析，在俯仰振荡中，飞机运动的数学描述为

$ \left\{ {\begin{array}{*{20}{l}} {\alpha (t) = {\alpha _0} + A{\rm{sin}}(\omega t) = {\alpha _0} + A{\rm{sin}}(2\pi ft)}\\ {q(t) = \dot \alpha (t) = 2\pi fA{\rm{cos}}(\omega 2\pi ft)}\\ {\bar q(t) = \frac{c}{{2V}}\dot \alpha (t) = {{\bar q}_{{\rm{max}}}}{\rm{cos}}(2\pi ft)} \end{array}} \right. $	（10）

式中：α₀为振荡的平均迎角；A为振幅；f为振荡频率；q为无因次俯仰角速率；q_max为一个周期振荡运动中最大无因次俯仰角速率，即

$ \left\{ {\begin{array}{*{20}{l}} {\bar q = \frac{{qc}}{{2V}}}\\ {{{\bar q}_{{\rm{max}}}} = \frac{c}{{2V}}A\omega = \frac{{\pi fAc}}{V}} \end{array}} \right. $	（11）

图 3给出了-4°≤α≤60°范围内不同振荡频率f、振幅A下俯仰阻尼导数试验结果，试验参数如表 2所示。

在常规飞行条件下，利用缩比模型开展动导数试验时，主要考虑减缩频率相似^[23]，以表征全尺寸飞机和模型之间关于刚体运动和绕流流体运动在时域上的同时性，形如：

$ \hat \omega = \left\{ {\begin{array}{*{20}{l}} {\frac{{\omega c}}{{2V}} = \frac{{\pi {f_C}}}{V}}&{{\rm{对于纵向运动}}}\\ {\frac{{\omega b}}{{2V}} = \frac{{\pi fb}}{V}}&{{\rm{对于横航向运动}}} \end{array}} \right. $	（12）

可见，常规试验中定义的减缩频率$\hat \omega $这个相似参数不包含关于运动振幅的信息，而从图 3给出的俯仰振荡试验结果看，在α≤10°飞机翼面流动未分离或α≥35°飞机翼面分离流动充分发展(或可认为流动已经完全分离)条件下，俯仰阻尼导数随减缩频率或振幅均不敏感，因此可以认为，在常规小迎角范围内或者在极限飞行状态中迎角特别大时，采用小振幅振荡获取动导数试验方法可以获得比较准确的动态气动力特性。但是，对于极限飞行状态中10° < α < 35°范围内(该范围通常为飞机失速到出现稳态尾旋之前)，俯仰阻尼导数不仅受减缩频率影响，即使飞机的减缩频率相同，运动的振幅也会对阻尼导数产生显著影响。可见，常规飞行条件下仅考虑减缩频率作为相似参数的试验方法对于部分极限飞行状态而言是难以捕捉其气动特征的。

图 4中分别给出小迎角(α₀=4°)、中大迎角(α₀=16°)、大迎角(α₀=40°)下，以振荡频率f=1 Hz、振幅A=5°开展动导数试验时，绘制的俯仰力矩系数C_m关于迎角的迟滞回线，图中不同“截止频率”指的是数据处理时的滤波截止频率。例如，截止频率为1 Hz(与强迫振荡频率相同，即通常说的保留1阶量)，即认为绕流流场变化频率与飞机刚体运动频率一致，迟滞回线为标准椭圆。众所周知，基于小扰动线性化假设的动导数就是由该1阶量计算得到的。可以看出，在常规小迎角范围或者在迎角特别大时，保留到高阶量(如6阶)与只留1阶量相比迟滞回线是比较接近的，即表示此时用动导数就可以比较准确地描述飞机的动态气动力，俯仰力矩随俯仰角速度是接近线性变化的；但是，如图 4(b)所示，在该中大迎角范围时，气动力存在明显的高阶量，因此，俯仰力矩随俯仰角速度变化有明显的非线性特征，此时，常规的线化导数不能精确反映动态气动力特性。从对气动力的频谱分析也可以看出，对应的数据含有1~6阶量(7阶以上基本没有)，一方面说明了该区域动态气动力存在显著非线性，另一方面该结果也表明对于气动力迟滞回线的计算，保留到6阶(当强迫振荡频率为1 Hz时，截止频率为6 Hz)是比较合理的。

上述结果说明了在10° < α < 35°范围内，式(7)所示的线性化假设不成立，必须对动态气动力的试验方法和数学模型进行改进。为了进一步分析动态气动力与角速率之间的关系，图 4(d)给出了振幅和频率均不同，但是在强迫振荡平均迎角处(α₀=15°)最大俯仰角速率相同的两条俯仰力矩迟滞回线。可见，尽管强迫振荡振幅和频率不同，但在最大角速率相同的位置，迟滞回线是非常接近的。这也是3.2节中把最大角速率引入到非线性动态气动力建模中的一个依据。

3.2 非线性气动力模型

由于在飞机的振荡运动中，角速率不仅与频率相关，而且振幅相关，式(11)定义的最大无因次角速率作为反映飞机动态运动的特征参数，可以同时体现频率和振幅对飞机动态气动特性的影响。为此，在本研究中引入最大无因次角速率作为阻尼导数的一个影响因子，以俯仰阻尼导数为例，即

$ {C_{m\bar q}}(t) = {C_{m\bar q}}(\alpha (t),\bar q(t)) $	（13）

用α₀=15°、A=10°的俯仰大振幅振荡为例对上述气动力模型进行分析和验证。俯仰大振幅运动如式(10)所示，根据表 2的动导数试验结果进行关于角速率的线性或非线性建模，线性数学模型表达式为

$ \begin{array}{l} {C_i} = {C_{i.{\rm{st}}}}(\alpha ) + {C_{i\bar q}}(\alpha )\frac{{qc}}{{2V}} + {C_{i\dot \alpha }}(\alpha )\frac{{\dot \alpha c}}{{2V}} = \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {C_{i.{\rm{st}}}}(\alpha ) + ({C_{i\bar q}}(\alpha ) + {C_{i\dot \alpha }}(\alpha ))\frac{{qc}}{{2V}} \end{array} $

（14）

非线性数学模型表达式为

$ \begin{array}{l} {C_i} = {C_{i.{\rm{st}}}}(\alpha ) + {C_{i\bar q}}(\alpha ,q)\frac{{qc}}{{2V}} + {C_{i\dot \alpha }}(\alpha ,q)\frac{{\dot \alpha c}}{{2V}} = \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {C_{i.{\rm{st}}}}(\alpha ) + ({C_{i\bar q}}(\alpha ,q) + {C_{i\dot \alpha }}(\alpha ,q))\frac{{qc}}{{2V}} \end{array} $

（15）

图 5、图 6给出了不同大振幅振荡频率(对应于不同的角速率区间)下，俯仰大振幅振荡试验结果和线性、非线性建模计算结果的对比。为了清楚识别图中气动力曲线特征，先对图中曲线形态、标注及其物理意义进行说明：以图 5为例，图中静态试验对应飞机俯仰角速度q=0(°)/s；在大振幅振荡运动中，由于动态气动力贡献，出现典型的“迟滞回线”，箭头表示对应曲线中飞机角运动的方向(箭头指向α增大的方向则q为正，反之为负)；对于俯仰力矩回线，气动力曲线随飞机运动沿着逆时针方向变化表明俯仰阻尼导数为负，飞机是动稳定的，反之则动不稳定，显然若回线出现交叉表明飞机动稳定性的突变，预示飞机可能会出现俯仰极限环振荡、机翼摇滚等非线性飞行力学现象。

可以看出，在飞机运动角速率相对较小区间内(图 5，|q|∈[0, 0.001 57])，线性模型和非线性模型计算结果基本一致，说明这种情况下飞机气动力非线性特征不明显。从两种模型计算结果和大振幅试验结果曲线对比来看，升力和俯仰力矩曲线拓扑结构一致，且俯仰力矩回线交叉点吻合，说明在该飞行状态参数区间内，式(14)和式(15)两个气动力模型能够捕捉飞机动稳定性关键特征：存在动稳定性突变点，在12° < α < 18°区域俯仰动不稳定。但与大振幅试验结果相比，迟滞回线幅值并不完全吻合。值得注意的是，上述两个气动力模型中，气动力大小都仅与当前时刻飞行状态参数(α和q)相关，而与其变化的时间历程无关，因而未蕴含气动力的非定常特征。

当角速率区间扩大(如图 6所示，|q|∈[0, 0.003 92])，气动力随俯仰角速率非线性变化特征凸显，线性气动力模型和非线性气动力模型计算结果有明显差异：若仍采用线性模型描述，俯仰力矩曲线拓扑结构仍保持与小角速率时一致，而大振幅试验结果显示，小角速度时的俯仰动不稳定区域已经退化到临界动稳定(f=1 Hz)或者动稳定状态(f=1.25 Hz)，常规的线化模型已经无法反映该特征；从非线性气动力模型曲线拓扑结构看，本文建立的非线性模型，则准确捕获了该动稳定性的非线性演化特征。同时，该情况下气动力非定常特性更为明显，进一步进行非定常气动力建模是必要的，这将在第4节进一步探讨。

综上所述，一方面，非线性模型可以更好捕获关于动稳定性的关键特征，因此，可同时结合减缩频率和最大无因次角速率(q_max)进行动导数试验设计和建模；另一方面，在量值上仍需进一步考虑非定常效应的影响。对于横航向进行分析也可得到类似的结果，本文中不再赘述。

4 动态气动力非定常特性

从第3节大振幅试验结果可以看出，加入非线性影响因素后，可以捕捉到飞机主要的动稳定性特征，但仅考虑非线性特性仍不足以精确定量描述运输机在进入极限飞行状态后可能面临的大迎角、快速运动时的气动特性，特别是在飞行失控运动中角速率较大时，气动力将呈现出较强的时间相关性。因此，在具体分析特定快速运动过程中的动态气动力时，还需考虑非定常特性影响。为了描述气动力的非定常特性，这里采用Goman-Khrabrov状态空间建模方法(G-K模型)^[10]，通过引入描述流场状态的内部变量即气流分离点位置，来描述分离流流场的动态发展过程，从而建立气动力响应的非定常模型。该模型已经成功应用在战斗机布局飞机非定常气动力建模中^{[15, 24]}。

本文以G-K模型以及Fan和Lutze等^[25]对该模型的发展为基础，以升力系数为例展开分析，模型表达式为

$ \left\{ {\begin{array}{*{20}{l}} {{\tau _1}\frac{{{\rm{d}}x}}{{{\rm{d}}t}} + x = {x_0}(\alpha - {\tau _2}\dot \alpha )}\\ {{C_L} = {C_{L0}} + {C_{L\alpha }}(x)\alpha + {C_{L\bar q}}(x)\bar q + {\Delta ^2}{C_L}} \end{array}} \right. $	（16）

式中：

$ {{x_0}(\alpha ) = \frac{1}{{1 + {{\rm{e}}^{\delta (\alpha - {\alpha _m})}}}}} $

$ {{C_{L\alpha }}(x) = {a_1} + {b_1}x + {c_1}{x^2}} $

$ {{C_{L\bar q}}(x) = {a_2} + {b_2}x + {c_2}{x^2}} $

$ {\Delta ^2}{C_L} = \frac{1}{2}[{C_{L{\alpha ^2}}}(x){\alpha ^2} + {C_{L{q^2}}}\left( x \right){\bar q^2} + 2{C_{L\alpha q}}(x)\alpha \bar q] $

其中：x∈[0, 1]表征流动分离点位置(对于全机布局飞机而言，该变量不再对应一个物理上的实际气流分离点，而认为是一个关于气流分离位置的等效描述)，x=1表示完全附着流，x=0表示完全分离流；时间常数τ₁表征的是流场改变后平衡分离流形成的过渡过程快慢，在气动特性上表现为动态迟滞回线跟踪静态曲线的趋势；时间常数τ₂表征的是流动分离和再附的迟滞，该迟滞与$\dot a$呈比例关系，在气动特性上反映为迟滞回线幅值大小；x₀(α)为定常状态下分离点的位置，α_m为分离点在中间时的迎角，该参数与δ一起确定静态时分离点位置随迎角的变化。可见，式(16)所示的非定常模型与式(15)的非线性模型相比，主要差别在于这里的气动导数依赖于内部状态变量x，从而依赖于流场结构及其变化的时间历程。

对于该非定常气动力模型，表征非定常特性的时间常数τ₁、τ₂的确定是关键。针对该问题，文献[10]推导了通过动导数来辨识时间常数的关系式，从而可以利用动导数试验数据确定时间常数。本文给出了另一种方法，即利用静态测力和大振幅振荡试验结果确定模型时间常数的方法，如下所述。

在式(16)的非定常模型中，只保留有关α、q的一阶项，忽略Δ²C_L及更高阶项，待确定的参数为：τ₁、τ₂、δ、α_m、C_L0、a₁、b₁、c₁、a₂、b₂、c₂。本研究分两步对参数进行辨识：

1) δ、α_m、C_L0、a₁、b₁、c₁只与静态气动力有关，因此，利用静态测力试验结果，用最小二乘拟合进行上述参数辨识。

2) 利用大振幅振荡试验结果，对τ₁、τ₂、a₂、b₂、c₂同样采用最小二乘拟合进行参数辨识。

作为一个计算实例，选择CRM飞机的一项升力静态测力试验结果和一项α₀=15°，A=10°，f=1 Hz的大振幅振荡试验结果，按上述方法进行参数辨识，结果为：τ₁=8.91 $\frac{c}{V}$，τ₂=2.53 $\frac{c}{V}$，δ=0.23，α_m=16.19，C_L0=0.05，a₁=2.67，b₁=0.23，c₁=2.54，a₂=-32.53，b₂=159.65，c₂=-139.25。

按照这组参数建立非定常气动力数学模型，与静态测力和大振幅试验结果对比如图 7、图 8所示。值得注意的是，上述参数只用了f=1 Hz 时的大振幅试验数据进行模型中与非定常特性有关的参数辨识，图 7给出了辨识的非定常气动力模型计算结果与试验结果对比，图 8中给出了使用该同一组参数(即τ₁、τ₂等参数保持不变)其他大振幅运动频率时的非定常气动力模型预测结果，与试验结果均吻合良好。可见，对该民机布局飞机，使用上述非定常气动力模型可以准确描述气动力的非定常特征。

综上，本研究表明，一方面，采用Goman-Khrabrov状态空间建模方法，可以准确描述典型民机布局飞机极限飞行状态下动态气动力的非定常特性，表明该模型适用于运输机布局飞机的非定常气动力建模，同时该模型基于流动机理进行建模，物理意义清晰，从而揭示了民机极限飞行状态下非定常气动力背后的流动结构演化和流动机理；另一方面，该研究为非定常模型时间常数的辨识提供了一种相对简便的通过大振幅振荡试验数据进行确定的方法。总体而言，利用第3节所述的非线性试验设计与建模方法可以获得飞机的主要动态气动力特征，预示出现非线性气动力的迎角范围，在此基础上，再针对特定的飞机运动过程，进一步开展大振幅或多自由度动态试验，进行非定常建模，可更准确获得运输机特定极限飞行状态的动态气动特性。该流程与方法为开展极限飞行状态动态气动力试验设计与建模研究提供了可行的技术途径。

5 结论

1) 在飞机超出常规迎角的极限飞行状态下，特别是在飞机开始失速到初始尾旋阶段(如对CRM布局而言，10° < α < 35°区域)，常规动态气动力试验的线性化、定常化假设不成立，需要考虑动态气动力的非线性与非定常特征。

2) 在动导数试验设计中，除了减缩频率，将一个振荡周期中的最大无因次角速率作为影响因素，形成关于迎角和角速率的非线性动态气动力模型，可以更好地捕捉飞机关于动稳定性的关键特征。

3) 采用Goman-Khrabrov状态空间建模方法，可以准确描述典型民机布局飞机极限飞行状态下动态气动力的非定常特性，表明该模型适用于运输机布局飞机的非定常建模。同时也揭示了民机极限飞行状态下非定常气动力背后的流动结构演化和流动机理。利用大振幅试验数据进行模型参数辨识，可以得到一种相对简单的通过动态试验确定Goman-Khrabrov模型中时间常数的方法。

4) 对于极限飞行状态的动态气动力研究，尽管动导数作为线性气动力参数已经不能用来描述飞机的运动形态，但仍可以预示非线性气动力的迎角范围。在此基础上，再针对特定的飞机运动过程，有针对性地进行大振幅试验及非定常气动力建模，可以更准确捕捉特定极限飞行状态的动态气动特性。因此，所建立的试验平台、方法、数据处理与建模等为民机极限飞行状态动态气动力风洞试验设计与建模研究提供了一个可行的途径。