跨声速非定常空气动力计算与分析

跨声速非定常空气动力计算Computation on Transonic Unsteady Aerodynamics北京大学力学与工程科学系理论与应用力学专业 00级陈雪梅摘要颤振问题一直是高速飞行器设计中的一大难题，特别在跨声速区段。

本文利用FLUENT6.1对一模型机翼的颤振行为进行了数值模拟，仿真机翼在高速气流中受激后扭曲变形最后发展成颤振的全过程，并对这一计算结果进行了初步分析，所得的算法具有普遍意义。

关键词：颤振，空气动力学，动网格[引言]早期的飞行器设计中的空气动力学分析都是将机翼﹑机身和其他气动部件当作刚体来处理。

但自第一架飞机诞生以来，空气动力学与飞机结构弹性的相互作用问题已经对航空技术的发展产生了重大影响，特别在‘彗星号’失事以后，人们对此倍加关心。

飞机在空气载荷作用下会出现可观的变形，这种变形将改变空气动载荷的分布，而它反过来又使变形发生变化。

在这种相互作用过程中，会引起振动，学术界称之为颤振。

这是一种自激振荡，它不断从气流中吸收能量。

当飞机发生颤振时，轻则出现不稳定和振动现象，重则因它引起材料‘疲劳’从而导致飞机在空中解体，以至机毁人亡。

在莱特兄弟首次试飞前，兰利的“空中旅行者”作了两次不成功的飞行试验。

第二次试飞时机翼和尾翼毁坏了，失败原因众说纷纭，气动弹性可能是第二次失败的罪魁祸首。

第一次世界大战中，英国的DH-9飞机尾翼颤振导致了飞行员死亡。

对此，英国空气动力学家贝尔斯托（L. Bairstow）首先对颤振进行了理论研究。

随着飞机速度的提高，空气动力增大，而重量小的结构形式使机翼抵抗变形的能力下降，所以气动弹性问题便得严重起来。

20世纪30年代初英国一家飞机连续发生有气动弹性引起的颤振事故，促使航空工程界对气动弹性问题普遍重视起来[摘自参考文献3，P118]。

其间的理论研究颇有成效。

美国力学家西奥多森（T. Theodorson）提交的研究报告对美国航空工业界建立颤振分析方法起了巨大作用。

单级跨音压气机三维非定常粘性流动计算西北工业大学刘前智周新海摘要本文以数值求解雷诺平均非定常Ｎ＿ｓ方程方法完成了单级跨音速压气机三维非定常粘性流动计算，给出了高分辨率ＨＨＤ格式的时间推进有限差分解．对某跨音单级压气机进行了定常和非定常流的详细数值计算，设计转速下定常粘性流场计算所得到的性能参数与实验结果吻合鞍好．在此基础上，对最高效翠状态进行了动／静叶相互子扰的非定常流动计算．非定常流计算模式由于更接近于真实流动，在高性能压气机的设计与分析中具有广阔的应用前景．关键词：非定常流动Ｎｓ：Ｙ程跨音速压气机动静叶相干一、引言在过去的十多年中所发展的三维定常流动的数值解法及其分析系统，已形成了目前的设计体系的基本框架，在工程实际中得到广泛应用．为了进一步地改进压气机和涡轮的气动设计方法，必须更准确地模述叶轮机中的实际流动．叶轮机中的流动是高度复杂的和非定常的，叶片不断承受着瞬时的、周期性的气动负荷及热负荷，这些非定常效应极大地影响叶轮机的气动性能和寿命．有效地预测叶轮机中动静叶排相干的非定常影响是改进目前的定常流设计方法和体系的一个重要课题，特别是对高性能多级跨音速风扇、压气机和涡轮来说，非定常干扰问题更加突出，追切需要发展相应的数值分析手段．国外从八十年代就开始了叶轮机动静叶排相干非定常流动方面的研究，从单级进展到研究多级叶排间的相干流动【１—４］．不过，早期的工作多集中在对涡轮流动的分析上，近来对风扇压气机的研究也逐渐增多．国内在这一领域的研究也已开始，主要集中于三维非定常无粘流动的分析研究上［５～６】．本文着重研究了跨音速轴流叶轮机动静叶排相干的三维非定常粘性流：＇ＪＪｌｉ３汁算问题，以无波动、无自由参数的耗散（ＮＮＤ）格式【７】为基础，得到了高分辨率格式非定常Ｎ．Ｓ方程时间推进有限差分解，保证了非定常流场中激波的分辨率．采用Ｌｕ．ＳＧＳ隐式解法［８Ｊ快速获得定常流场解，以最大限度地减少计算时间．利用本文方法对某单级跨音速压气机进行了数值分析，首先进行了设计转速下的不同状态定常粘性流场计算．在此基础上，选取了最高效率状态进行了非定常流动计算，并对计算结果进行了分析．本文计算所得到的周期性非定常解，有效地预估了动静叶排相干非定常效应及其对叶轮机内部流动的影响．二、基本方程、边界条件将圆柱坐标系（ｐ，ｒ，ｚ）下的三维非定常粘性流动的平均Ｎ·Ｓ方程转换到一般曲线坐标系（ｆ，玎，ｆ）下，其形式为固／ａ＋反Ｅ—Ｅｖ）／ａ；＋反Ｆ一易）／却＋文Ｇ—Ｇ，）／Ｇ＝Ｈ（１）式中各项的意义详见文献（９］．分子粘性系数ｕ，由Ｓｕｔｈ口＇ｌａｎｄ公式计算，而紊流粘性系数ｕｔ，由Ｂａｌｄｗｉｎ－Ｌｏｍａｘ模型【ｌｏ】确定．进口边界（转子进口）上，给定总压、总温、气流的绝对切向和径向纾，＊度沿径向的分布，出口边界（静子出口）上，给定内径处的静压，出口其他径向位置处的压力利用径向平衡方程确定．在物面边界，为减少网格数，采用壁面函数［１１】求物面剪切应力，这时，在物面应用了速度滑移条件．周期边界的处理：当动静叶数目不相等（甚至相互不成整数倍）时，选择ｍ个转子叶８７片通道与ｎ个静子叶片通道，ｍ和ｎ韵选取原则是尽量减少静子叶片几何调整的程度，又不致使计算工作量成倍增加．此时，转子叶排与静子叶排将用组合栅距（含数个叶片通道）来取代原单通道的栅距，在一个组合栅距上流动参数满足周期性条件．动静叶排问信息的传递将主要取决于叶排间交界面的处理方法．在定常流计算中，采用周向平均模型．在非定常流动计算时，则采用直接插值的方法传递叶排间的信息．三、数值方法．本文采用如下方法进行方程的时间离散：在作为非定常计算的初始场计算中，采用ＬＵ－ＳＧＳ隐式解法，得到定常解．为了保证时间精度，非定常流动计算则采用了显式方法，其时间步长满足稳定性条件．方程的空间离散采用高分辨率ＮＮＤ格式。

低超声速来流中机翼跨声速非定常气动力数值解法
张建柏;张秀滨
【期刊名称】《空气动力学学报》
【年(卷),期】1993(011)003
【总页数】5页(P303-307)
【作者】张建柏;张秀滨
【作者单位】不详;不详
【正文语种】中文
【中图分类】V211.41
【相关文献】
1.跨声速机翼非定常气动力的全位势粘位迭代计算 [J], 陆志良;Voss.,R
2.超声速机翼—尾翼组合体的非定常气动力边界元法 [J], 周文伯;裘松刚
3.超声速来流中横向喷流角度对流动与混合特性的影响 [J], 杨银军;窦志国;段立伟
4.带操纵面机翼的非定常跨声速流数值解法 [J], 余涛; 张建柏
5.跨声速机翼操纵面非定常气动力Euler方程计算 [J], 董璐;郭同庆
因版权原因，仅展示原文概要，查看原文内容请购买。

超声速翼型气动力对比计算研究超声速翼型是指在超过声速的空气流动条件下，流经翼型表面的速度大于声速。

超声速翼型的气动特性与亚声速翼型有很大的不同，因此对超声速翼型的气动力进行研究和计算，对于超声速飞行器的设计和性能优化具有重要意义。

超声速流动下的气动力分为横向力、升力和阻力三个方向。

首先，横向力是指垂直于飞行方向的力，用于控制超声速飞行器的横侧稳定性。

当超声速飞行器在飞行过程中发生偏离时，通过调整横向力的大小和方向可以使其重新回到预定的飞行轨迹。

横向力的计算可以通过数值模拟方法，如CFD（计算流体力学）进行。

通过CFD计算，可以得到超声速流动下的横向力系数。

其次，升力是指垂直于飞行方向且指向上方的力，用于支持超声速飞行器的重量。

超声速翼型的升力计算一般采用低阻力曲线翼型理论方法或CFD方法。

低阻力曲线翼型理论方法是根据翼型的几何形状和气动力学性质，通过一系列的计算公式推导得出升力系数。

CFD方法则是通过数值模拟方法对超声速流动的速度场和压力场进行计算，然后根据经验公式和实验数据计算升力系数。

最后，阻力是指沿飞行方向的阻碍运动的力，它是超声速飞行器在飞行中需要克服的主要力量。

超声速流动下的阻力包括压力阻力和摩擦阻力。

压力阻力是由于超声速流动中气体与翼型表面的冲击作用以及翼型附近湍流引起的。

摩擦阻力是由于超声速流动中气体与翼型表面的摩擦力引起的。

计算超声速流动下的阻力需要考虑这两种阻力的贡献，可以通过CFD方法进行计算，也可以采用经验公式进行估算。

综上所述，超声速翼型气动力的计算研究对于超声速飞行器的设计和优化至关重要。

通过数值模拟方法如CFD，可以计算得到超声速流动下的横向力系数、升力系数和阻力系数，为超声速飞行器的性能分析和改进提供重要依据。

同时，可以通过低阻力曲线翼型理论方法和经验公式估算得到升力系数和阻力系数，为初期设计提供快速的气动力估算手段。

基于CFD气动力辨识模型的气动弹性数值计算聂雪媛;刘中玉;杨国伟【摘要】Unsteady aerodynamic reduced order model (ROM)based on computational fluid dynamics (CFD)can either improve computational efficiency or hold the same computational accuracy as CFD.However,the modeling of ROMs based on the identification technology is limited by the loading form and/or frequency spectrum of excitation signals.To overcome the shortcoming and improve the identification efficiency,the random white noise signal was taken as an excitation to model unsteady aerodynamic forces with one-mode-at-a-time loading form.A multiple-input multiple-output (MIMO)model was obtained by the linear superposition of identified single-input multiple-output(SIMO)autoregressive moving average (ARMA)models.With excitation signals of different frequencies and shape on the identified model,the simulation results by ROM and direct CFD computation were compared and it indicates that the same computation accuracy as by CFD can be acquired by ROM.Coupled with structural model,the ROM has been used to predict the flutter boundaries of the wing of AGARD445 .6 .The numerical simulations show that the flutter results predicted by ROM are in general agreement with those by direct computation of unsteady Navier-Stokes equations and wind-tunnel experiments.It verifies that ROM can provide a high efficient method for transonic aeroelastic numerical analysis.%基于系统辨识技术的非定常气动力降阶模型（Reduced-OrderModels，ROMs），保留流体力学（Computa-tional Fluid Dynamics，CFD）计算精度的同时能提高计算效率，但对输入信号加载方式或/和信号频谱要求较高。

速度场演化结构非定常气动力速度场演化结构非定常气动力是指空气流动的速度场在时间和空间上都发生了改变，造成气动力的变化。

在实际的气动力学问题中，流体的速度场通常是非定常的，因为空气在流动过程中，会受到外界环境的干扰和影响，以及流体本身的特性，例如流体粘性和湍流等。

因此，研究速度场演化结构非定常气动力对于设计和优化飞行器、建筑物以及其他需要考虑气动力的工程问题具有重要意义。

1. 运动方程：空气流动的速度场可以通过Navier-Stokes方程来描述。

Navier-Stokes方程是流体力学中最基本的方程之一，它描述了流体的连续性、动量守恒和能量守恒。

在非定常条件下，Navier-Stokes方程需要引入时间项来描述速度场的演化过程。

2. 湍流模型：湍流是非定常流动的一种常见形式，它在风力发电、飞行器设计等领域具有重要应用价值。

湍流的处理是气动力学研究中的一个重要问题，需要使用湍流模型来描述湍流的统计性质。

常用的湍流模型有雷诺平均Navier-Stokes方程（RANS）和大涡模拟（LES）等。

3.边界条件：边界条件对于非定常气动力的研究至关重要。

边界条件包括速度场、压力场以及物体表面的边界条件。

在非定常情况下，需要考虑流动的时间变化和空间变化对边界条件的影响。

对于速度场演化结构非定常气动力的研究，可以采用数值模拟的方法进行。

数值模拟是一种通过利用计算机对气动力学问题进行数值计算和分析的方法。

数值模拟方法中常用的有有限元法、有限体积法和有限差分法等。

在实际应用方面，速度场演化结构非定常气动力的研究对于飞行器的设计和安全性评价具有重要意义。

例如，对于飞机的气动外形设计和燃油效率的优化，需要考虑机翼和机身等部件对速度场的影响以及速度场演化对气动力的变化。

此外，速度场演化结构非定常气动力的研究还可以应用于建筑物的抗风设计、风洞实验以及地质工程等领域。

总而言之，速度场演化结构非定常气动力是一个复杂的问题，涉及到流体力学、湍流模型和数值模拟等多个领域的知识。

第十六届全国计算流体力学会议跨声速机翼突然失速现象的 CFD 研究吴佳莉，马海，卜忱（中航工业空气动力研究院，哈尔滨市第 88 号信箱 150001) 摘 要：20 世纪 90 年代中后期，F/A-18E/F 在跨声速机动飞行中多次出现了非指令横向运动引起的“掉翼尖”现象，使飞机的横向稳定性和操纵性下降，对飞行安全构成极大的威胁。

本文基于中航工业空气动力研究院自主研 发的 UNSMB 非结构网格计算平台对带有前缘锯齿的某先进布局飞机基本构型开展了跨声速机翼突然失速 （AWS） 特性计算研究， 分析了引起 AWS 现象的复杂流动机理， 达到了初步预测机翼突然失速特性的能力。

研究结果显示，1阶段原型机进行试飞时，技术人员发现该飞机在马赫数 0.7~0.95，高度 3000~12000 米，迎角 6°~12°的跨 声速飞行条件下，一侧的机翼会发生突然失速（AWS） ，诱发另一侧机翼下沉，伴随非指令性的横向运动， 称之为“掉翼尖”现象，飞机的横向稳定性和操纵性下降，对飞行安全构成极大的威胁[1]。

这个问题后来由 波音公司相继通过改变前缘襟翼偏转程序和在机翼折叠机构上增开了一个多孔门得到初步解决。

AWS 问题 一度使 F/A-18E/F 的研制试飞工作陷于停滞，因此，在 20 世纪 90 年代末美国空军启动了“机翼突然失速” （AWS）专项研究项目，以求借助 CFD 计算、 风洞试验、 地面模拟和飞行试验等多种手段， 了解 F/A-18E/F 和其他飞机布局 AWS 现象发生的机理，并发展用于评估和预测机翼突然失速的计算品质因素。

AWS 现象的实验研究主要是通过对跨声速风洞静态/动态试验得到的力、力矩、压力和 PSP 图像进行 研究得到。

McMillin[2]等人在 LaRC 16ft 跨音速风洞对 8%缩比的预生产型 F/A-18E 模型的 AWS 现象进行了 系统研究（见图 1） 。

风力机风轮非定常气动载荷计算1.引言随着清洁能源技术的发展，风能作为一种可再生、清洁的能源被广泛应用于电力生产领域。

而风力机作为转化风能为电能的设备，其稳定性和可靠性对于电力系统的稳定运行具有重要作用。

然而，风力机受到非定常风速和风向的影响，导致风轮非定常载荷，影响其稳定性和可靠性。

因此，研究非定常气动载荷计算方法对于风力机运行的控制和优化具有重要意义。

2. 非定常气动载荷特点2.1 风力机非定常风场特点风力机非定常载荷来源于风场的非定常性和风轮本身的非定常性。

其中，风场的非定常性是由于风速和风向的变化导致的，而风轮本身的非定常性则是由于风轮运动状态的变化引起的。

风速变化包括风向变化、风速周期性变化、突然风暴等。

这些变化导致风力机受到的非定常载荷具有以下特点：（1）涡旋生成：当风速和风向发生变化时，会在风轮背风侧产生涡旋，引起非定常载荷变化。

（2）波动载荷：风速周期性变化会引起非定常载荷的周期性变化。

（3）外加载荷：风暴风等突然变化的风速和风向变化会引起较大的外加载荷。

2.2 风力机非定常气动载荷特点风力机非定常气动载荷是指风轮运动状态变化引起的载荷变化。

风轮运动状态的变化包括旋转角速度的变化、叶片变形等。

而这些变化会导致风轮的气动载荷发生变化，具有以下特点：（1）非定常气动力：当风轮旋转时，气动力也随着变化。

这种气动力具有特殊的非定常特性，例如相位滞后、自激振荡等。

（2）非定常扭矩：风轮非定常气动力的变化会引起扭矩的变化，这种非定常扭矩会对风力机的稳定性和可靠性产生影响。

（3）振动载荷：风轮非定常气动载荷的变化会引起风轮的振动，这种振动载荷会对风力机的结构强度和寿命产生影响。

3. 非定常气动载荷计算方法为了控制和优化风力机的运行，需要对其受到的非定常气动载荷进行计算和分析。

目前，非定常气动载荷的计算方法包括解析方法、半经验方法和试验方法。

3.1 解析方法解析方法是一种基于物理原理和数学模型的计算方法，可以计算出理论上的非定常气动载荷。

飞机T型尾翼跨音速颤振特性研究杨飞;杨智春【摘要】由于飞机T型尾翼的结构与气动布局特点,T型尾翼颤振计算不能套用常规尾翼的分析方法,而需要考虑平尾面内运动以及静升力等因素的影响.而跨音速空气压缩性效应和非定常气动力计算的不准确性,使得T型尾翼跨音速颤振计算更加困难,准确性较低.因此,需要采用试验为主计算为辅的方法来研究飞机T型尾翼跨音速颤振特性.针对某T型尾翼结构,用ZAERO软件等价片条势流跨音速颤振(ZTAIC)方法计算T型尾翼跨音速颤振特性,研究了马赫数、风洞气流密度和平尾迎角对T 型尾翼颤振特性的影响.通过升力系数斜率空气压缩性修正计算方法和跨音速颤振模型风洞试验方法得到了飞机T型尾翼的跨音速颤振的凹坑曲线和空气压缩性特性,两种方法得到结果一致.【期刊名称】《振动与冲击》【年(卷),期】2013(032)010【总页数】5页(P50-54)【关键词】跨音速;颤振;T型尾翼;风洞;试验;跨音速凹坑;压缩性【作者】杨飞;杨智春【作者单位】中国商飞上海飞机设计研究院强度部,上海200232;西北工业大学航空学院,西安710072【正文语种】中文【中图分类】V215.3飞机跨音速颤振特性从根本上决定了飞机的颤振包线，事关飞机稳定性安全。

通常可以通过试验或计算手段得到飞机的颤振临界耦合模态、临界颤振动压、跨音速颤振动压压缩性系数和颤振裕度。

在亚音速(低马赫数)情况下，空气压缩性对颤振速度影响较小，当马赫数大于0.5时，必须考虑空气压缩性的影响，在马赫数等于1.0附近的跨音速区，颤振速度(颤振动压)会急剧降低，形成一个所谓“跨音速凹坑”。

飞机T型尾翼是指平尾位于垂尾稍部，平尾和垂尾组成一个“T字”结构形式的尾翼。

T型尾翼结构具有诸多优点，一方面，T型尾翼布局可使平尾避开机翼尾流或尾吊发动机喷流的影响，增大平尾力臂、提高操纵效率;另一方面，T型尾翼构型可以实现后机身大开口，便于大型装备的货物装运，同时T型尾翼的高置平尾可满足水上飞机设计要求。

气动力学研究中的非定常问题分析与计算方法研究气动力学是研究物体在气流中运动时所产生的各种力和运动规律的学科。

研究对象包括飞行器、导弹、火箭、汽车等。

气动力学的发展在军事、航空、航天、汽车等领域具有重要的应用价值。

在气动力学研究中，非定常问题分析与计算方法研究是重要的研究方向。

一、非定常问题的定义非定常问题是指在气体流动过程中，流场参数随时间变化的问题。

对于这类问题，需要采用有关的计算方法和分析技术。

二、非定常问题的研究现状非定常问题的研究是气动力学领域的热点之一，现有研究成果较多，包括计算方法和分析技术等方面。

研究成果广泛应用于实际项目，提高了相关领域的技术水平。

1、计算方法非定常问题的计算方法主要有数值方法和解析方法两种。

数值方法主要包括有限差分法、有限元法、谱方法、蒙特卡罗方法等。

解析方法主要是基于数学推算，采用解析技巧求解非定常问题。

这些计算方法在不同领域的气动力学问题中得到广泛应用。

2、分析技术非定常问题的分析技术主要包括实验研究和理论分析两种方法。

实验方法主要是通过模型试验和现场试验来获取相关数据，并分析数据来研究非定常问题。

理论分析方法主要是通过推理和数学分析方法研究非定常问题的规律和特性。

三、计算方法的主要研究内容非定常问题的计算方法是研究非定常问题的重要手段。

近年来，国内外学者对计算方法进行了深入研究，主要内容包括以下方面：1、非定常流动的模拟和仿真计算方法是非定常问题研究的重要手段，其中，非定常流动的模拟和仿真是目前的研究热点。

国内外学者通过数值计算和实验研究等手段，对非定常流动规律进行了深入探究，取得了一些重要成果。

2、非定常流动的数值模拟在非定常问题的数值模拟研究中，数值方法是实现非定常流动模拟的重要手段。

近年来，学者们对数值方法进行了深入研究，如有限体积法、有限元法、伴随法等方法等，这些方法在实践中得到了广泛应用，为非定常问题研究提供了重要支撑。

3、非定常问题的策略和技术非定常问题是一个复杂的问题，需要采用多种策略和技术实现研究。

第41卷第11期2020年11月哈㊀尔㊀滨㊀工㊀程㊀大㊀学㊀学㊀报Journal of Harbin Engineering UniversityVol.41ɴ.11Nov.2020NACA0012翼型跨声速强迫运动非定常气动力模型张庆1,2,叶正寅3(1.西安航空学院飞行器学院,陕西西安710077;2.南洋理工大学机械与航空工程学院,新加坡639798;3.西北工业大学航空学院,陕西西安710072)摘㊀要:传统的一阶线性叠加的气动力模型不再适用于现代高机动性飞行器的非定常气动力建模,为了考察更高阶的气动力模型对非定常迟滞效应模拟的适用程度,本文采用自主发展的求解器,分别计算了NACA0012翼型在跨声速来流条件下做单自由度强迫沉浮㊁俯仰以及沉浮/俯仰两自由度耦合运动的非定常气动力的变化规律㊂然后在Etkin 气动力模型的基础上,探讨了不同类型的高阶的气动导数在非定常气动力建模中的作用㊂研究结果表明:将Etkin 气动力模型中升力和俯仰力矩对迎角的导数项由一阶拓展至二阶就可以较为精确地重构出翼型在强迫运动各阶段的非定常升力和俯仰力矩㊂关键词:跨声速;气动导数;气动力建模;沉浮运动;俯仰运动;耦合运动;Etkin 模型;非定常气动力DOI :10.11990/jheu.201903018网络出版地址:http :// /kcms /detail /23.1390.u.20201028.1424.010.html 中图分类号:V211.41㊀文献标志码:A㊀文章编号:1006-7043(2020)11-1683-06Unsteady aerodynamic model of NACA 0012associated with forcedoscillations and translations in transonic flightZHANG Qing 1,2,YE Zhengyin 3(1.School of Aircraft,Xiᶄan Aeronautical University,Xiᶄan 710077,China;2.School of Mechanical &Aerospace Engineering,Nan-yang Technological University,Singapore 639798,Republic of Singapore;3.School of Aeronautics,Northwestern Polytechnical Uni-versity,Xiᶄan 710072,China)Abstract :It is unsuitable to model unsteady aerodynamics for high-agility modern aircraft by the traditional first-or-der linear superposition theory.For this study,in order to investigate the applicability of a higher order aerodynam-ic model to the simulation of hysteresis effects,the unsteady time histories of aerodynamics for NACA0012associat-ed with single-freedom forced motions,plunging and pitching,and coupled plunging and pitching motion under transonic conditions were investigated computationally based on in-house codes.The effects of various aerodynamic derivatives on aerodynamic models are discussed based on the Etkin aerodynamic model.Final results indicated that unsteady lift and pitching moment in forced single or coupled motions could be accurately regenerated if the Etkin model is expanded to the second order derivative of the angle of attack with respect to time.Keywords :transonic;aerodynamic derivative;aerodynamic model;plunging;pitching;coupled motion;Etkin model;unsteady aerodynamics收稿日期:2019-03-06.网络出版日期:2020-10-28.基金项目:国家自然科学基金重点项目(11732013);校级科研基金项目(2018KY1226).作者简介:张庆,男,讲师,博士;叶正寅,男,教授,博士生导师.通信作者:张庆,E-mail:zhangqing2220@.㊀㊀气动导数作为描述飞行器机动飞行和受扰动时气动特性的关键性气动参数,在飞行器气动性能㊁控制系统和总体设计中扮演着非常重要的作用[1-4]㊂在传统的飞行动力学相关问题的研究中,气动力的数据往往基于小扰动线性叠加原理计算出来,在这种准定常假设情况下,气动力仅仅表示为瞬时飞行状态参数的函数,并且可以以一种简单的解析函数关系式表示出来[2-5]㊂但是,现代飞行器的飞行包线普遍向大迎角区域扩展,在大迎角下飞机机动飞行产生的三维非定常分离流和涡流使得空气动力呈现高度非线性特性,气动力和力矩不仅依赖于瞬时迎角㊁侧滑角㊁姿态角等参数,而且与它们的时间历程有关,因此原来使用的低阶线性叠加模型将不再适用[5-6]㊂同时,由于机动飞行状态涵盖了较大的迎角㊁侧滑角㊁角速率的变化范围,如果采用风洞实验或是数值计算模拟,其时间成本和经济成本都难以接受[7-10]㊂哈㊀尔㊀滨㊀工㊀程㊀大㊀学㊀学㊀报第41卷因此有必要建立起较大飞行包线内普适性较好的的非定常气动力模型[1,4]㊂Etkin模型是目前动导数求解时最常用的一种非定常气动力模型,Etkin模型物理意义明确,考虑了时间历史效应对气动导数的影响[3]㊂但是,在非定常气动力建模时,该模型中的各项气动导数对不同运动形式的非定常气动力的影响规律和适用程度尚不清楚㊂为此,本文结合Etkin气动力模型,研究了气动力关于迎角的一阶和二阶导数在气动力模型的作用,希望能精确地重构出翼型单自由度或是耦合强迫运动过程中的非定常气动力,为未来发展高效的㊁可靠的气动力模型提供参考数据㊂1㊀强迫运动非定常气动力模型本文的计算采用课题组自己开发的柔性体动力学问题求解软件GMFlow[11-13],其中流场求解部分采用基于SA模型的有限体积法[13],强迫运动时的网格变形方法为弹簧网格变形方法[14-16]㊂为了验证求解方法的正确性,首先计算了标准算例NACA0012翼型强迫俯仰运动的非定常气动力变化情况,将计算结果与文献中的计算结果和实验结果对比,对比结果见文献[13]㊂俯仰运动的运动规律可以描述为[15]:α(t)=α0+A sin(ωt)=α0+A sin(2πft)(1)式中:α0是初始位置处的迎角;A是简谐振动的振幅;ω是简谐振动的圆频率;f是简谐振动的频率㊂本文定义减缩频率为:k=ωC2Vɕ(2)式中C是翼型的弦长㊂在本文中,强迫运动时自由来流的马赫数为0.755,翼型弦长为1.0m,强迫运动的减缩频率为0.0814㊂俯仰运动的初始迎角为0.016ʎ,俯仰振幅为2.51ʎ㊂图1(a)是强迫俯仰运动时的升力系数和关于1/4弦点的俯仰力矩系数随时间的变化曲线,图中计算了3个周期的气动力,由图可知,在第1个计算周期的初始阶段,计算的结果收敛性较差,这主要是由于定常计算的步数不足㊂从第2个周期开始,力和力矩系数已经达到了较好的谐振性,可以认为计算结果已经收敛㊂因此,为了减小计算量,本文的所有强迫运动过程都只计算了3个运动周期㊂图1(b)是翼型强迫沉浮运动时的力和力矩系数变化情况,其运动规律为:z(t)=z0+z m sin(ωt)(3)式中:z0=0是初始位置处的纵向位移;z m=0.1m 是沉浮运动的振幅㊂考虑洗流影响,在沉浮运动的任一时刻,瞬时迎角为:α(t)=α0-ωz m cos(ωt)/Vɕ(4)㊀㊀图1(c)是翼型强迫俯仰/沉浮耦合运动时的力和力矩系数变化情况,其运动规律为式(1)和式(3)叠加㊂对比图1可知,虽然耦合运动形式是俯仰和沉浮运动的叠加,但是耦合运动的气动力和力矩并不等于俯仰运动和沉浮运动的简单叠加,这也说明了翼型强迫运动时气动力的非线性迟滞特性比较复杂,并不是简单的线性叠加关系㊂图1㊀不同运动过程升力和力矩系数随时间变化Fig.1㊀History of lift/moment coefficients in different mo-tions㊃4861㊃第11期张庆,等:NACA0012翼型跨声速强迫运动非定常气动力模型1.1㊀一阶气动模型根据Etkin 气动力模型[2-3],强迫运动过程中的非定常气动力可以表示为:ΔC j =C j -C j 0=C jαΔα+C j ̇αC 2V ɕ()Δ̇α+C jq C 2V ɕ()Δq (5)式中C j 0是平衡位置处的力系数或是力矩系数㊂由于式(5)中C j ̇α和C jq 的量纲相同,都是空气动力系数对角度随时间一阶变化率的导数,所以在本文中称式(5)为一阶气动力模型㊂根据强迫俯仰运动时运动规律可知:̇α(t )=q =Aωcos(ωt )(6)㊀㊀俯仰运动的非定常气动力可以表示为:ΔC j =A sin(ωt )㊃C jα+kA (cos(ωt )-1)(C j ̇α+C jq )(7)㊀㊀此处需要注意,由于ΔC j 是相对于初始位置的变化量,因此右侧是(cos(ωt )-1)而不是cos(ωt )㊂所以:C j ̇α+C jq =ʏT n +1T nΔC j cos(ωt )d t()kAπω()(8)㊀㊀根据强迫沉浮运动时运动规律可知:̇α(t )=ω2z m sin(ωt )/V ɕ(9)㊀㊀沉浮运动的非定常气动力可以表示为:ΔC j =C jα-ωz m cos(ωt )V ɕ+C j ̇αω2z m V ɕC2V ɕ()sin(ωt )(10)㊀㊀所以:C j ̇α=ʏT n +1T nΔC jsin(ωt )d t ()ω2z m πC2V 2ɕ()(11)㊀㊀将式(8)和式(11)相减就可以得到单独的C j ̇α和C jq ,根据图1(a)和图1(b)的结果,可以求得一阶气动力模型的各个导数值,升力系数对̇α和q 的导数值分别为-38.9457和6.6747,力矩系数对̇α和q 的导数值分别为-2.0595和-1.3344㊂1.2㊀二阶气动模型根据Etkin 气动力模型[2-3],非定常气动力可以表示为:ΔC j =C jαΔα+C j ̇αD 2V ɕ()Δ̇α+C jα㊃㊃D 2V ɕ()2Δα㊃㊃+C jqD 2V ɕ()Δq +C j ̇qD 2Vɕ()2Δ̇q (12)㊀㊀由于式中C jα㊃㊃和C j ̇q 的量纲相同,都是空气动力系数对角度随时间二阶变化率的导数,所以在本文中称式(12)为二阶气动力模型㊂与1.1节类似,由俯仰运动的气动力变化规律可以得到:C j ̇q=C jα/k 2-ʏT n +1T n ΔC j sin(ωt )d t ()k 2A πω()C j ̇α+C jq =ʏT n +1TnΔC j cos(ωt )d t ()kA πω()ìîíïïïïïï(13)㊀㊀由沉浮运动的气动力变化规律可以得到:C j ̇α=ʏT n +1T nΔC jsin(ωt )d t ()ω2z m πC 2V ɕ2()C jα㊃㊃=ʏT n +1T nΔC j cos(ωt )d t()(z m πV ɕ()+C jα)/k 2ìîíïïïïïï(14)㊀㊀式(13)减去式(14)就可以得到单独的C j ̇α和C jq ,根据图1(a)和图1(b)的结果,可以求得一阶和二阶气动力模型的各个气动导数值,一阶导数与上节完全相同,升力系数对α㊃㊃和̇q的导数值分别为578.5118和-38.2752,力矩系数对α㊃㊃和̇q的导数值分别为18.8996和12.9044㊂式(10)㊁(11)与式(12)~(14)相比,一阶气动导数完全一样,这也间接说明传统上忽略高阶导数的做法对低阶气动导数的求解结果并没有影响,这是传统上普遍采用Etkin 气动力模型进行小迎角㊁小扰动飞行包线范围内动态稳定性分析的重要原因㊂2㊀气动力建模结果比较为了定量考察这些气动力模型对强迫运动过程非定常迟滞效应模拟的适用程度,本节对比了这些气动力模型的计算结果与直接采用CFD 进行计算得到的结果㊂图2分别是采用一阶和二阶Etkin 气动力模型计算得到的强迫俯仰运动㊁强迫沉浮运动以及耦合运动的气动力与采用CFD 方法得到的气动力迟滞曲线的对比图㊂由图2(a)可知,对于强迫俯仰运动,采用二阶气动导数得到的升力系数与CFD 计算值完全重合,而采用一阶气动导数得到的升力系数误差随着迎角的增加而增大,在最大迎角位置比CFD 计算值大50%㊂对于俯仰力矩系数,一阶模型的误差较大,而二阶模型的结果与CFD 计算值虽然不像升力系数那样完全重合,但是吻合程度也较好㊂表1是俯仰运动不同位置处的不同变量对该时刻非定常气动力的贡献情况,需要注意的是,强迫俯仰运动时̇α与q 的数值相等,α㊃㊃与̇q的数值相等㊂由表1可知,在俯仰运动1/4周期时,到达抬头最大位置处,此时α㊃㊃对应的非定常升力和力矩分别占总的非定常升力和力矩的-37.98%和-45.72%,̇q对应的非定常升力和力矩分别占总的非定常升力和力矩㊃5861㊃哈㊀尔㊀滨㊀工㊀程㊀大㊀学㊀学㊀报第41卷的2.51%和-31.21%㊂由于忽略了α㊃㊃与̇q的影响,所以一阶气动力模型对应的升力和力矩系数与CFD 计算值差别较大,见图2(a)㊂在俯仰运动3/4周期时,到达低头最大位置处㊂此时α㊃㊃与̇q这2项对非定常特性的贡献也比较大,所以与1/4周期时的情况类似,一阶模型的计算结果误差较大㊂而在1/2周期时,翼型低头经过初始位置,由于Δα㊁Δα㊃㊃以及̇q的数值刚好为0,此时非定常气动力主要由̇α和q 项产生,所以一阶气动力模型就能较为准确地重现出非定常气动力㊂在一个周期时,翼型抬头经过初始位置,此时的非定常气动力贡献情况与1/2周期时相反,此时由于Δ̇α和Δq 的数值刚好为0,此时非定常气动力主要由α㊁α㊃㊃以及̇q项产生㊂这些分析结果与图2(a)的结论一致,说明在俯仰运动过程中,二阶模型才能更准确地再现出非定常气动力和力矩㊂图2㊀不同运动过程升力和力矩系数迟滞曲线Fig.2㊀Comparison of lift /moment coefficients indifferent motions表1㊀俯仰运动不同位置非定常气动力分布情况Table 1㊀Percentage distributions at different time in thepitching motion %周期α̇αα㊃㊃q ̇q 0.25T 0.50T 0.75T 1.00TΔC L 109.4431.41-37.98-5.38 2.51ΔC m 76.0861.20-45.7239.65-31.21ΔC L0120.680-20.680ΔC m 060.68039.320ΔC L 228.27-65.52-79.2211.23 5.24ΔC m -74.8160.1844.9538.9930.69ΔC L 147.950-51.350 3.40ΔC m -8935.8505369.583666.27㊀㊀图2(b)是采用一阶和二阶Etkin 气动力模型计算得到的强迫沉浮运动的气动力与采用CFD 方法得到的气动力迟滞曲线的对比图㊂由图2(b)可知,对于强迫沉浮运动,采用二阶气动导数得到的升力系数与CFD 计算值几乎重合,而采用一阶气动导数得到的升力系数误差较大,在最大纵向位移位置比CFD 计算值大90%㊂对于俯仰力矩系数,一阶模型的误差较大,而二阶模型的结果与CFD 计算值虽然不像升力系数那样完全重合,但是吻合程度也较好㊂表2是沉浮运动不同位置处迎角的各阶导数对非定常气动力的贡献情况,由于强迫沉浮运动时没有俯仰角速度,所以在表2中没有出现q 和̇q 项对应的非定常气动力㊂由表2可知,在沉浮运动1/4周期时,到达纵向最大位置处,此时α㊃㊃对应的非定常升力和力矩分别占总的非定常升力和力矩的-94.84%和148.24%㊂由于忽略了α㊃㊃的影响,所以一阶气动力模型对应的升力和力矩系数与CFD 计算值差别较大,见图2(b)㊂在沉浮运动3/4周期时,到达纵向最小位置处㊂此时α㊃㊃对非定常特性的贡献也比较大,其数值分别为-36.92%和-49.93%,所以与1/4周期时的情况类似,一阶模型的计算结果误差较大㊂而在1/2周期时,翼型下沉经过初始位置,由于Δ̇α的数值刚好为0,此时非定常气动力主要由α和α㊃㊃产生㊂在一个周期时,翼型回到初始位置,此时的非定常气动力贡献情况与1/2周期时刚好相反,由于Δα和Δα㊃㊃的数值为0,所以此时非定常气动力主要由Δ̇α产生㊂这与图2(b)的结论一致,说明在沉浮运动过程中,二阶模型才能更准确地再现出非定常气动力和力矩㊂图2(c)是采用一阶和二阶Etkin 气动力模型计算得到的强迫俯仰/沉浮耦合运动的气动力与采用CFD 方法得到的气动力迟滞曲线的对比图㊂由图2(c)可知,对于强迫耦合运动,采用二阶气动导数得到的升力系数与CFD 计算值完全重合,而采用一阶气动导数得到的升力系数误差较大,在最大纵向位移位置比CFD 计算值大145%㊂对于俯仰力矩系数,一阶模型的误差较大,而二阶模型的结果与CFD 计算值虽然不像升力系数那样完全重合,但是㊃6861㊃第11期张庆,等:NACA0012翼型跨声速强迫运动非定常气动力模型吻合程度也较好㊂表2㊀沉浮运动不同位置非定常气动力分布情况Table 2㊀Percentage distributions at different time in theplunging motion %周期α̇αα㊃㊃0.25T 0.50T 0.75T 1.00TΔC L273.27-78.43-94.84ΔC m -246.70198.46148.24ΔC L153.150-53.15ΔC m 250.570-150.57ΔC L 106.3930.53-36.92ΔC m 83.0966.84-49.93ΔC L 0100.000ΔC m 0100.000㊀㊀表3是耦合运动不同位置处各导数对非定常气动力的贡献情况㊂由表3可知,在耦合运动1/4周期时,到达纵向最大位置处,此时也处于抬头的最大位置,此时α㊃㊃对应的非定常升力和力矩分别占总的非定常升力和力矩的-45.35%和-70.82%,̇q 对应的非定常升力和力矩分别占总的非定常升力和力矩的2.19%和-35.25%㊂由于忽略了α㊃㊃与̇q 这2项的影响,所以一阶气动力模型对应的升力和力矩系数与CFD 计算值差别较大,见图2(c)㊂在耦合运动3/4周期时,到达纵向最小位置处,此时也处于低头的最大位置㊂此时α㊃㊃对非定常特性的贡献也比较大,其数值分别为-245.66%和21.26%,̇q 对非定常特性的贡献分别为25.87%和23.00%,所以与1/4周期时的情况类似,一阶模型的计算结果误差较大㊂而在1/2周期时,翼型下沉经过初始位置,由于Δ̇q的数值刚好为0,对非定常气动力没有贡献㊂在一个周期时,翼型回到初始位置,此时Δq 的数值刚好为0,对非定常气动力没有贡献㊂虽然耦合运动气动力并不是单独运动的简单叠加,但是通过对经典Etkin 气动力模型的二阶延拓,能准确地再现出沉浮/俯仰耦合运动过程的非定常气动特性㊂表3㊀耦合运动不同位置非定常气动力分布情况Table 3㊀Percentage distributions at different time in thecoupled pitching /plunging motion %周期α̇αα㊃㊃q̇q0.25T 0.50T 0.75T 1.00TΔC L 130.6617.18-45.35-4.68 2.19ΔC m 117.8643.43-70.8244.78-35.25ΔC L77.3459.74-26.84-10.240ΔC m 25.2154.58-15.1535.360ΔC L 707.86-443.48-245.6655.4125.87ΔC m -35.2261.8521.1629.2123.00ΔC L 172.08-16.19-59.930 4.04ΔC m -349.6492.22210.81146.613㊀结论1)不论是强迫俯仰运动㊁沉浮运动,还是俯仰/沉浮耦合运动,将气动导数拓展至迎角和俯仰角的二阶导数,都可以十分精确地重现出强迫运动过程中的非定常升力变化情况㊂2)由于俯仰力矩的迟滞曲线并不是简单的椭圆形,二阶模型计算出的强迫运动过程的俯仰力矩与CFD 计算值的吻合程度不像升力那么好,说明俯仰力矩的模型要比升力更加复杂㊂3)俯仰/沉浮耦合运动的非定常气动力并不是俯仰运动和沉浮运动的简单叠加,说明精确的气动力建模还需要深入考虑其他变量的影响㊂本文的研究结果表明,Etkin 气动力模型对于非线性较强的气动力建模仍然具有较好的适用性,但是,对于三维流动以及接近失速迎角情况下的非定常气动力的建模,需要更加深入地讨论马赫数㊁减缩频率㊁更高阶导数以及交叉导数在非定常气动力模型中的作用㊂参考文献:[1]杨磊,叶正寅.倾转涵道倾转过渡阶段的非定常气动力[J].航空动力学报,2015,30(1):155-163.YANG Lei,YE Zhengyin.Unsteady aerodynamic force oftilt ducted fan during transition period[J].Journal of aero-space power,2015,30(1):155-163.[2]张庆,叶正寅.基于气动导数的类X -37B 飞行器纵向稳定性分析[J].北京航空航天大学学报,2020,46(1):77-85.ZHANG Qing,YE Zhengyin.Longitudinal stability analysis for X -37B like trans-atmospheric orbital test vehicle based on aerodynamic derivatives[J].Journal of Beijing Universi-ty of Aeronautics and Astronautics,2020,46(1):77-85.[3]ETKIN B.Dynamics of atmospheric flight[M].Mineola:Dover Publications,2012.[4]袁先旭,陈琦,谢昱飞,等.动导数数值预测中的相关问题[J].航空学报,2016,37(8):2385-2394.YUAN Xianxu,CHEN Qi,XIE Yufei,et al.Problems innumerical prediction of dynamic stability derivatives [J].Acta aeronautica et astronautica sinica,2016,37(8):2385-2394.[5]和争春,钱炜祺,朱国林,等.飞行器跨声速区俯仰力矩系数建模方法研究[J].空气动力学学报,2005,23(4):470-475.HE Zhengchun,QIAN Weiqi,ZHU Guolin,et al.Re-search on the modeling of pitching moment coefficient in transonic condition for flight vehicle[J].Acta aerodynami-ca sinica,2005,23(4):470-475.[6]GHOREYSHI M,CUMMINGS R M.Challenges in the aer-odynamics modeling of an oscillating and translating airfoil at large incidence angles[J].Aerospace science and tech-nology,2013,28(1):176-190.[7]NELSON R C,PELLETIER A.The unsteady aerodynamicsof slender wings and aircraft undergoing large amplitude㊃7861㊃哈㊀尔㊀滨㊀工㊀程㊀大㊀学㊀学㊀报第41卷maneuvers[J].Progress in aerospace sciences,2003,39 (2/3):185-248.[8]MARZOCCA P,LIBRESCU L,SILVA W A.Aeroelastic response of nonlinear wing sections using a functional series technique[J].AIAA journal,2002,40(5):813-824.[9]SILVA W A.Application of nonlinear systems theory to transonic unsteady aerodynamic responses[J].Journal of aircraft,1993,30(5):660-668.[10]SILVA W A.Simultaneous excitation of multiple-input/multiple-output CFD-based unsteady aerodynamic systems [J].Journal of aircraft,2008,45(4):1267-1274.[11]ZHANG Qing,HUA Ruhao,YE Zhengyin.Experimentaland computational investigation of novel vertical tail buffet suppression method for high sweep delta wing[J].Science China technological sciences,2015,58(1):147-157.[12]ZHANG Qing,YE Zhengyin.Novel method based on in-flatable bump for vertical tail buffeting suppression[J].Journal of aircraft,2015,52(1):367-371. [13]张庆.高速再入飞行器动力学问题研究[D].西安:西北工业大学,2018.ZHANG Qing.Research on flight dynamics of high-veloci-ty reentry vehicles[D].Xiᶄan:Northwestern Polytechni-cal University,2018.[14]RAUSCH R D,BATINA J T,YANG H T Y.Spatial ad-aptation of unstructured meshes for unsteady aerodynamic flow computations[J].AIAA journal,1992,30(5): 1243-1251.[15]BATINA J T.Unsteady Euler airfoil solutions usingunstructured dynamic meshes[J].AIAA journal,1990, 28(8):1381-1388.[16]KANDIL O A,CHUANG H A.Unsteady transonic airfoilcomputation using implicit Euler scheme on body-fixed grid[J].AIAA journal,1989,27(8):1031-1037.本文引用格式:张庆,叶正寅.NACA0012翼型跨声速强迫运动非定常气动力模型[J].哈尔滨工程大学学报,2020,41(11):1683-1688. ZHANG Qing,YE Zhengyin.Unsteady aerodynamic model of NACA0012associated with forced oscillations and translations in transonic flight[J].Jour-nal of Harbin Engineering University,2020,41(11):1683-1688.㊃8861㊃。

3.6.6 亚音速、跨音速和超音速飞行以及气动力系数的变化1亚音速、跨音速和超音速飞行图3-40所示为随着飞行马赫数的提高，机翼翼表面上激波变化的情况，从图中可以看出，当M a=0.72时，翼型上表面首次出现了等音速点，这个翼型的临界马赫数M a临=0.72。

当M a=0.77时，在翼型上表面首次出现了局部超音速区和局部激波，激波分离也可能在这时出现。

随着M a数继续提高，等音速点向前移，局部激波向后移。

超音速区逐渐扩大。

当M a=0.82时，下翼面开始出现局部激波。

随着Ma数的继续提高，翼型表面的超音速区继续扩大，直到M a=1.05 ，局部激波移到了翼型的后缘，在翼型的前缘形成了脱体正激波，这时，只有在正激波的后面有一块亚音速区，其他流场已全部变成超音速了。

如果继续提高Ma数，亚音速区会进一步缩小，大约在M a=1.3时，就可以认为气流在翼型表面全部都是超音速流动了。

图3-40 随着马赫数M a的增加，激波逐渐产生（1）亚音速飞行：在飞行M a《M a临（一般为0.7左右）时，气流流过机翼表面的流场全部都是亚音速流场，在这个范围内，飞机的飞行是亚音速飞行。

（2）跨音速飞行：在飞行M a>M a临,在机翼表面出现了局部超音速区和局部激波后，直到机翼流场全部称为超音速流场之前（M a临<M a≤1.3），这个范围内飞机的飞行是跨音速飞行。

飞机进行跨音速飞行时，机翼表面的流场既有亚音速流场又有超音速流场。

（3）超音速飞行：到飞行M a>1.3 以后，机翼表面的流场全部成为超音速流场，飞机的飞行就是超音速飞行了。

2.随着飞行M a数的提高，气动力系数的变化随着飞行飞行M a数的提高，翼型表面的流场发生着剧烈的变化，翼型的空气动力也着发生变化。

图2-41所示为升力系数C L、阻力系数C D以及交点X随着M a数提高而变化的情况。

从图中可以看到，从M a>M a临开始，位置FX约为25%左右，并基本保持不变。