第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㊃。
