基于CFD气动力辨识模型的气动弹性数值计算

聂雪媛，刘中玉，杨国伟

（中国科学院 力学研究所 流固耦合系统力学重点实验室，北京 100190）

对飞行器气动弹性不稳定性预测在设计中非常重要。随计算流体力学（Computational Fluid Dynamics，CFD）、计算结构力学（Computational Structure Dynamics，CSD）及计算机性能的发展，直接耦合CFD／CSD的数值模拟方法成为非线性气动弹性分析中可信度最高方法［1］。然而在气动弹性研究中，非定常气动力计算占据大部分计算时间，使基于CFD的气动弹性模拟难以在实际工程应用。如何既能保证计算精度又能提高计算效率成为气动弹性的研究热点。非定常流场模型降阶技术 （ReducedOrder Models，ROMs）的出现为CFD／CSD数值模拟方法用于工程设计成为可能。降阶模型由全流场全阶CFD模型近似投影获得低阶模型，该模型能捕捉非线性流场动力学特征，计算效率远高于CFD。非定常气动力ROMs代替CFD与结构模型耦合，能进行气动弹性分析或气动伺服弹性分析及设计［2－3］。气动力降阶模型基于动态线性假设［4］，在给定非线性定常流场基础上利用小扰动理论推导获得到线性模型。在基于系统辨识的CFD降阶模型方法中，Silva等［5］提出基于脉冲响应的Volterra级数模型降阶方法，并与CFL3D结合进行气动弹性颤振分析。但此方法需获得系统全部Markov参数，而所用辨识数据因通过“一次激励结构一个模态”获得，会造成结构模态个数较多时所需脉冲响应个数呈幂函数增长。为此，文献［6－7］开发出同时激励结构多个模态方法提高辨识效率。认为输入多个激励信号间须存在相位差；但未给出如何确定相位差，需多次调试。Gupta等［8－10］提出采用飞行试验中易实现的速度为3211信号作为输入进行系统辨识。利用自回归滑动平均模型（Autoregressive Moving Average，ARMA）作为替代CFD的降阶模型。该方法亦需对各模态同时输入的3211信号进行反复试验、调整频率，以期输出响应能含希望激发的频带。为进一步提高效率，Raveh［11］用多个白噪声信号模拟结构多模态同步输入，辨识三种多输入、输出（Multiple Input Multiple Output，MIMO）模型，即频域模型、ARMA模型及状态空间模型，分析AGARD445.6机翼在跨声速区域的颤振边界。该方法虽可避免对输入信号频谱范围要求，但仍未涉及多个输入同步信号间相位差的确定。徐敏等［12－13］以一次激励一个模态方式的阶跃响应信号为输入，基于Volterra级数方法建立非定常气动力降阶模型，并用于二维机翼气动弹性主动控制及气动弹性分析。张伟伟等［14］用3211信号建立基于CFD的 ARMA模型用于 Isogai机翼及AGARD445.6机翼气动弹性分析，并认为输入信号需反复调试。刘晓燕等［15］采用基于观测器的Volterra核辨识方法，通过每个模态单独激励气动力系统方式，在每个模态输入含结构感兴趣频率正弦信号建立气动力降阶模型用于气动弹性非线性气动力预测及颤振抑制。该方法本质仍属基于Volterra级数的ROMs。

为方便输入信号选择、避免同步输入多信号时反复测试方可确定相位差、提高辨识计算效率，本文以随机白噪声信号为系统输入，用一次激励一个模态方式同时运行多个CFD数值模拟程序计算系统响应。利用所得多组训练数据辨识出多个基于ARMA的单输入多输出（Single Input Multiple Output，SIMO）部分气动力降阶模型。因非定常气动力建模建立于动态线性假设理论基础上，故利用线性系统叠加原理获得基于CFD的多输入、输出非定常气动力降阶模型，并将其与结构模型耦合用于AGARD445.6机翼颤振边界预测。通过将预测结果与风洞试验结果比较表明，本文所提建模方法准确、有效。

1 基于ARMA的气动力降阶模型

颤振边界预测为气动弹性分析中的重要计算及研究弹性体在小扰动激励下系统响应，符合动态线性模型分类［16］，因此可用线性模型描述非定常气动力，进行气动弹性颤振分析。

本文所用气动力降阶模型进行气动弹性颤振边界分析过程可描述为将给定的结构模态位移一次一个模态输入到非定常流场求解器，利用N－S方程求出各阶模态对应的广义气动力，设关心结构的前n个模态，则会获得n组“单输入－n输出”（记为1－n）的模型训练数据，用系统辨识技术可得n个“1－n”的线性模型，利用线性系统叠加原理，最终获得“n个输入－n个输出”（记为n－n）的气动力降阶模型，见图1下方虚线框。用所得降阶模型替代CFD与结构模型耦合，在结构振幅小扰动下通过改变动压，可快速获得结构位移响应，从而迅速找到某马赫数下的颤振临界参数，见图1上方虚线框。此过程所耗时间与CFD／CSD耦合计算相比可忽略。

图1 基于气动力辨识模型的气动弹性分析Fig.1 Aeroelastic analysis based on aerodynamic force identification model

本文所用气动力降阶模型为自回归滑动平均ARMA模型，即时间域离散差分方程，形式为

式中：ya为n个广义气动力输出组成的列向量；u为单模态位移输入信号标量；na，nb为输出、输入延迟阶数；k为离散时间变量；A，B为欲辨识系统模型参数。

设结构模态数为n（以下n均为结构前n阶模态），则式（1）所描述的为图1下方虚线框中每个“单模态位移输入－n个模态广义气动力输出”的1－n ARMA模型。据文献［14］取4～8阶延迟模型即可精确描述非定常效应，因此利用MATLAB的arx函数可较快辨识最好的表征训练数据（式（1））最优的1－n ARMA模型。获得1－n ARMA模型后用文献［17］方法在Matlab中编程实现将离散差分方程转换为离散状态空间方程形式。第i阶模态激励下状态空间方程为

上述系数矩阵均为分块矩阵。式（2）为第i阶模态激励下所得“单输入－n输出”的SIMO系统状态空间方程。为便于气动弹性分析，需获得“n个模态输入激励下n个模态广义气动力输出”的多输入多输出MIMO系统状态空间模型。因此对n个 “单输入－n输出”状态空间方程线性叠加，获得所需MIMO状态空间方程为

各参数中上标为模态阶次。

2 激励信号

模型辨识主要为激励信号选取，要求所用激励信号能激起被辨识系统所有频段信号。随机信号为含频带范围最宽的信号，理论上用其作为激励信号能获得系统所有频段范围内响应。为实际气动弹性工程应用，仅研究某特定频率附近模态运动所致气动力，对随机信号进行带通滤波处理。本文主要研究AGARD445.6机翼前四阶模态 9.60 Hz、38.16 Hz、48.35 Hz、91.54 Hz引起的广义气动力。对输入信号进行滤波，保证激励信号包含该频段即可。

输入信号类型确定后，数据采样频率及长度将决定辨识精度。由于该激励信号用于CFD计算，故其采样时间与CFD计算时间步相同，本文为无量纲时间0.05，对应有量纲频率 fs为 11 523.6 Hz（马赫数为0.96时），其 Nyquist频率为5 761 Hz，远高于所需最高频率91.54 Hz。对白噪声激励信号采用带通方式滤波，保留频段范围为0.000 1fs～0.1fs。采样长度据频率分辨率定义需在CFD上运行约2 500个时间步，即可辨识系统最低频率9.60 Hz。因此可获得精度较高模型。虽运行2 500个时间步需耗费一定计算时间，但与频繁修改参数进行CFD计算及气动弹性分析的传统方法相比，计算效率改善仍相当可观。

图2 机翼近壁面网格图Fig.2 Mesh grid of near the wing

3 算例与分析

为验证该辨识方法的有效性及准确性，本文对AGARD445.6进行模型辨识，利用辨识所得模型进行颤振边界、颤振速度预测，并与风洞试验结果比较。计算网格单元数约40万。近壁面网格见图2。CFD网格计算中首层网格厚度约为根弦长的8×10－5，结构振型归一化见文献［18］。采用非定常 N－S方程计算机翼弹性运动时非定常流场，空间采用基于结构网格技术的有限体积法，选HLLE格式进行空间离散，用代数无限插值法实现网格运动。用双时间步进行时间离散，拟时间步为隐式格式LUSGS方法。施加模态位移激励时需保证产生的结构实际位移、速度满足线性扰动假设。

本文采用高斯白噪声信号，在输入CFD前先进行滤波处理。所选带通范围 0.000 1fs～0.1fs，其中 fs为CFD计算用采样频率。对应有量纲带通约1.152～1 152 Hz。理论上为辨识出最低频率，需至少2 500个时间步的CFD计算结果作为模型训练数据，但实际应用中的ARMA模型，用200个时间步足够，多余数据可用于校验辨识的模型。由于针对结构前四阶模态，因此只要辨识四个“单输入－4输出”的ARMA模型。模型辨识结果阶次na＝4，nb＝5。

图3为马赫数0.96、输入第一阶模态位移为过滤高斯白噪声信号时辨识模型与CFD计算所得气动力见图3，其中虚线为模型计算结果，实线为CFD计算结果。两种方法计算结果吻合。

图4为马赫数0.96、给定第二阶模态位移50 Hz正弦信号时辨识所得ARMA模型计算的气动力与CFD直接计算的气动力结果比较见图4，其中虚线为模型计算结果，实线为CFD计算结果，可见二者吻合较好。

图3 第一阶模态激励下辨识模型与CFD计算结果对比Fig.3 Comparison of ARMA model simulated responses and CFD response to excitation of the 1st mode

图4 第二阶模态激励下辨识模型与CFD计算结果比较Fig.4 Comparison of ARMA model simulated response and CFD response to sinusoidal excitation of the 2nd mode

图5 输入3211信号Fig.5 3211 Excitation of the first four modes

三、四阶模态激励的模型准确性检验与一、二阶类似。在检验辨识出的四个“单输入－4输出”ARMA模型准确性基础上对其分别进行状态空间转换获得（式（2））方程再进行线性叠加，最终获得“4输入－4输出”气动力状态空间模型。对该空间模型精确度采用3211信号进行检验。输入具有一定相位差的4个3211信号见图5。在图5信号输入下“4输入－4输出”状态空间模型计算的广义气动力与CFD直接计算结果比较见图6，其中实线为CFD计算结果，方块为“4输入－4输出”状态空间模型计算结果。由图6看出，MIMO系统状态空间模型计算结果与CFD直接计算结果吻合较好为评价辨识模型计算精度，本文用常用评价模型准确度参数即模型拟合度为

图6 模型输出与CFD计算结果比较Fig.6 Comparison of identification model simulated results and CFD results

式中：为CFD计算结果；为辨识模型计算结果；m为数据长度；η＝1为ARMA模型与N－S方程计算结果完全一致。

白噪声信号、50 Hz正弦信号及3211信号作为结构模态位移输入时，ARMA的模型拟合度见表1。由表1看出，所辨识模型精度较高。

表1 不同输入信号激励下模型拟合度Tab.1 Comparison of model fitness to different excitation signals

采用辨识模型预测的AGARD445.6机翼颤振频率边界及速度边界见图7、图8。由二图看出，耦合结构的N－S方程及ROM计算结果与风洞试验结果均较吻合，且机翼颤振速度在跨声速区域出现特有的“凹坑”现象。表明本文所用非定常气动力模型辨识方法能用于跨声速气动弹性稳定性分析。

图7 AGARD 445.6机翼无因次颤振频率边界Fig.7 Nondimensional flutter frequency boundary of AGARD 445.6 wing

图8 AGARD 445.6机翼无因次颤振速度边界Fig.8 Nondimensional flutter speed boundary of AGARD 445.6 wing

基于系统辨识模型进行气动弹性分析效率较高。如寻找马赫数为0.96时的颤振临界参数。采用求解非定常N－S方程方法搜寻一次参数至少需CFD运行2000个周期，若每个周期含20次内迭代（1次内迭代约3 s），运行时间约需2000 min（Intel Xeon 2.4 GHz双CPU），而改变参数重新计算仍费时2 000 min。用非定常气动力建模技术时以产生2 500个辨识用数据而言，需CFD运行2 500 min，而参数辨识时间在Matlab下约需10 s。故模型确定后搜寻一次参数，Matlab仅需1～2 s。

4 结 论

（1）本文以过滤随机白噪声信号为激励信号，利用动态线性假设及线性系统叠加原理将辨识所得MISO系统ARMA模型线性叠加获得MIMO系统模型，实现非定常气动力建模。将辨识所得模型用白噪声信号、正弦信号、3211信号等进行检验，并与直接CFD计算结果比较，从而验证辨识模型的准确性。该方法与现有辨识方法相比，具有对激励信号要求少、易于施加、辨识效率高等优点。

（2）将辨识所得时域信号变换为状态空间模型用于AGARD445.6机翼颤振边界预测。气动力降阶模型与结构模型耦合计算结果与CFD／CSD耦合计算结果及风洞试结果一致，表明本文方法不仅效率显著提高，能用于气动弹性稳定性分析，并为气动伺服弹性分析及设计提供可能。

［1］陈刚，李跃明.非定常流场降阶模型及其应用研究进展与展望［J］.力学进展，2011，41（6）：686－701.CHEN Gang，LI Yueming.Advances and prospects of the reduced order model for unsteady flow and its application［J］.Advances in Mechanics，2011，41（6）：686－701.

［2］Dowell E H.Eigenmode analysis in unsteady aerodynamics：Reducedorder models［J］.AIAA Journal，1996，34（8）：1578－1583.

［3］Lucia D J，Beran P S，Silva W A.Reducedorder modeling：new approaches for computational physics［J］.Progress in Aerospace Sciences，2004，40（1／2）：51－117.

［4］Dowell E H.A modern course in aeroelasticity，3rd Revised and Enlarged Edition［M］.New York：Klewer Academic Publishers，1995.

［5］Silva W A.Bartels R E.Development of reducedorder models for aeroelastic analysis and flutter prediction using the CFL3Dv6.0 code［J］.Journal of Fluids and Structures，2004，19：729－745.

［6］Silva W A. Simultaneous excitation of multipleinput multipleoutput CFDbased unsteady aerodynamic systems［J］.Journal of Aircraft，2008，45（4）：1267－1274.

［7］Kim T，Hong M，Bhatia K G，et al.Aeroelastic model reduction for affordable computational fluid dynamicsbased flutter analysis［J］.AIAA Journal，2005，43：2487－2495.

［8］Gupta K K，Voelker L S，Bach C，et al.CFDbased aeroelastic analysis of the X43 hypersonic flight vehicle［C］.Proceedings of the 39th Aerospace Sciences Meeting and Exhibit，2001.

［9］Cowan T J，Gupta K K.Accelerating CFDbased aeroelastic predictions using systems identification［C］.36th AIAA Atmoshperic Flight Mechanics Conference and Exhibit，1998：85－93.

［10］Cowan T J，Gupta K K. Development of discretetime aerodynamic model for CFDbased aeroelastic analysis［C］.Proceedings of the 37th Aerospace Sciences Meeting and Exhibit，1999.

［11］Raveh D E. Identification of computationalfluiddynamic based unsteady aerodynamic models for aeroelastic analysis［J］.Journal of Aircraft，2004（41）：620－632.

［12］徐敏，陈刚，陈士橹.基于非定常气动力低阶模型的气动弹性主动控制律设计［J］.西北工业大学学报，2004，22（6）：748－752.XU Min，CHEN Gang，CHEN Shilu.Design method of active control law based on reducedorder model of unste ady aerodynamics for aeroelasticity［J］.Journal of Northwestern Polytechnical University，2004，22（6）：748－752.

［13］姚伟刚，徐敏.基于Volterra级数降阶模型的气动弹性分析［J］.宇航学报，2008，29（6）：1711－1716.YAO Weigang，XU Min.Aeroelastic numerical analysis via Volterra series approach［J］.Journal of Astronautics，2008，29（6）：1711－1716.

［14］张伟伟，叶正寅.基于非定常气动力辨识技术的气动弹性数值模拟［J］.航空学报，2006，27（4）：579－583.ZHANG Weiwei，YE Zhengyin.Numerical simulation of aeroelasiticity basing on identification technology of unsteady aerodynamic loads［J］.Journal of Astronautics，2006，27（4）：579－583.

［15］刘晓燕.基于非定常气动力降阶模型的气动弹性研究［D］.北京：北京航空航天大学，2011：54－61.

［16］Dowell E H，Crawley E F.A modern course in aeroelasticity［M］.Kluwer Academic，1989.

［17］Cowan T J， Gupta K K. Development of discretetime aerodynamic model for CFDbased aeroelastic analysis［R］.AIAA，1999.

［18］Carson Y E.Agard standard aeroelastic configurations for dynamic responseiwing 445.6［C］.AGARD－R－756，1988.