大振幅振荡来流条件下非定常气动力模型计算验证与弱可压缩性修正

郭 力, 吕计男, 冯 峰, 王 强

(中国航天空气动力技术研究院, 北京 100074)

大振幅振荡来流条件下非定常气动力模型计算验证与弱可压缩性修正

郭 力, 吕计男, 冯 峰, 王 强*

(中国航天空气动力技术研究院, 北京 100074)

针对大振幅振荡来流条件下薄翼受到的非定常气动力，Isaccs和Greenberg分别发展了非定常气动力模型。这两种模型可以用于直升飞机桨叶与风力发电叶片的气动力分析，模型在不可压缩无黏来流条件下建立，但实际流动不可避免粘性和弱可压缩性的影响，需要检验两种模型的适用性。针对粘性效应的影响，2014年Strangfeld对于NACA0018翼型，通过风洞实验验证了在Reynolds数25万时，Isaccs和Greenberg的模型仍适用，实验的Mach数为0.0326，流动近似不可压缩流动。针对可压缩性的影响，通过数值模拟方法进行了研究。首先重复了实验在Mach数为0.0326时的结果，并进一步考察了当Mach数提高为0.1、0.2和0.3时非定常气动力的变化。结果表明随着Mach数的提高，升力系数的最高点逐渐高于模型，并且相位逐渐落后，在Mach数为0.3时差别最明显，非定常升力系数最高点计算与模型相差50%。此即表明弱可压缩性对模型的预测结果影响不可忽略。为了扩展模型在Mach数变化时的适用范围，对模型进行了弱可压缩性修正。通过考虑速度变化引起均匀来流中密度的变化，修正了翼型附近流体密度，使其跟随来流Mach数变化。采用此方法，将计算与模型的幅值差别减小到5%左右。

Isaacs模型；Greenberg模型；非定常气动力；低马赫数；弱可压缩性修正

0 引 言

旋转的叶片, 如直升飞机的旋翼, 风力机的叶片等, 即使在定常运动时, 由于旋转的作用, 其截面受到的来流大小与迎角也会非定常变化。针对翼型的平动与非定常迎角变化问题, Theodorsen在1935年提出了非定常气动力的理论解[1]。Isaacs根据Theodorsen理论加入了来流速度变化的影响[2]与迎角变化的影响[3]。Greenberg考虑了速度脉动与翼型运动的复合作用[4]。Mateescu在2006年发展了在低速可压缩情况下非定常气动力模型[5], 模型采用在翼型前缘与脊(ridge)上的速度奇点的解, 考察了刚性与柔性翼型的平动与转动非定常性。Walker等考虑了翼型在柔性变形情况下的气动力[6]。Ramesh发展了模型用于考察前缘涡对于昆虫飞行的影响[7]。Cho等采用3维非定常涡格(vortex lattice)方法, 发展了非定常气动力模型, 用于考虑速度脉动引起的非定常气动力对气动弹性稳定性的影响[8]。刘启宽等[9]采用Greenberg公式分析了风力机叶片的挥舞特性。Williams等采用Greenberg模型[4]设计控制率, 抑制来流振荡引起的受力变化[10]。Jose等采用Theodorsen的理论分析了二维多段翼型间隙对非定常气动力的影响[11]。Strangfeld[12]使用NACA0018翼型通过实验对Iscaccs和Greenberg提出的模型进行了考察。实验中翼型固定, 来流速度做正弦变化。来流的平均速度为11m/s, 最大速度为16.5m/s, 最小速度为5.5m/s。在海平面条件下, 平均速度对应马赫数近似为0.0326, 流动可以近似为不可压缩流动。基于弦长的Reynolds数为2.5×105。实验考察了迎角为2°与8°时升力系数随速度的变化。在迎角为2°时实验结果振动较大。但通过前两阶Fourier模态的拟合结果与理论吻合较好。在迎角为8°时实验较光滑, 但低估了最大升力系数。这是由于迎角8°时实验中涡在翼型表面的分布与理论分布不同引起。

Isaacs和Greenberg的气动力公式采用势流理论推导。其试用范围为: 二维不可压缩势流流动, 迎角α较小, 满足近似关系α≈sin(α), 来流不会反向。本文通过计算, 考察弱可压缩性对Isaacs和Greenberg模型精度的影响。文中首先重复了Strangfeld等的实验条件下升力系数随速度的变化。然后对更高马赫数0.1, 0.2和0.3时升力系数的变化进行考察。

1 Greenberg、Isaacs非定常升力系数气动力模型

来流速度围绕平均速度做正弦振动时会引起翼型受到的力和力矩也随着速度变化, 做正弦变化。设来流速度u做正弦变化的表达式为:

Cl、Cm和lm的表达式为无穷级数,Strangfeld等[12]将Cl中的无穷级数在第8项之后(m>8)截断;lm的级数表达式在第15项(n>15)后截断。其误差小于10-6。

根据Greenberg模型, 非定常升力系数Cl(t)与定常升力系数Cl,st之比为:

式中a是取矩的点距离前缘的距离。在不可压缩无粘势流假设下, 对于翼型四分之一弦长的地方取矩, 无论迎角多少, 得到的力矩系数近似为0[16]。文中所有力矩均关于翼型的前缘取矩。

图1为Iscaacs和Greenberg模型的比较。其中实线为Isaacs模型结果, 虚线为Greenberg模型结果。无论对于升力系数, 还是力矩系数,Greenberg公式在最低点与最高点处低估了升力。但是在初始位置与接近平衡点时两个模型结果相近。

2 翼型绕流计算设置

Isaacs模型与Greenberg模型均对于不可压缩流动提出。为了与模型的适用范围接近, 本文计算采用的马赫数较低。分别考察了马赫数为0.0326、0.1、0.2和0.3时, 迎角为2°和8°的情况。其中马赫数为0.0326近似为Strangfeld实验[12]的来流条件。计算采用HUNS3D程序[13]。通量部分采用AUSM+-up格式[14], AUSM+-up格式在马赫数较低时可以保证格式的收敛性和精度[14]。湍流模型采用SA模型[15], Reynolds数为Re=2.5×105。时间推进采用双时间步方法, 其中隐式迭代采用LU-SGS方法。计算采用2维网格, 计算区域为矩形, 如图2所示。

在流向与法向方向上为40倍的翼型弦长, 出口的上下两角导圆。翼型位于计算区域中央。迎角通过翼型旋转设定。来流方向始终与入口垂直, 上下两个边界采用滑移无穿透边界条件, 使得入口下游的流场能够根据来流速度的变化进行调整, 避免人为指定带来的误差。入口边界指定速度, 其随时间的变化为式(1)。式(1)中σ=0.5与Strangfeld等的实验[12]所用振幅相同,ω=4.72×10-3根据Strangfeld的实验[12]中约化频率k=0.0074得到。

为了验证网格与计算条件的适用性, 采用程序对来流马赫数为0.0326、迎角为8°的定常情况进行计算, 并与实验进行对比。图3为翼型表面升力系数。其中实线为计算所得结果, 点为实验结果。在翼型上表面与下表面上计算与实验吻合较好。在前缘处, 下表面得到的升力系数与实验接近。在前缘上表面和尾缘处压力略高于实验所得值, 但差别在误差范围内。

3 非定常气动力特性分析

3.1 模型与计算比较

为了验证程序在非定常时的计算结果, 将计算得到的升力系数Cl、力矩系数Cm在马赫数0.0326时与实验和理论进行了比较。为了方便比较, 升力系数采用定常状态做了归一化。非定常计算从定常状态启动, 为了消除启动时初始效应对计算的影响, 结果为速度变化第二个周期时的数据。图4、图5中实线为Isaacs的理论结果, 虚线为Greenberg的理论结果，点划线为计算结果, 实心圆点为Strangfeld等的实验结果[12]。实验中压力传感器受到噪声的影响, 出现脉动[12]。为了减小实验时压力传感器噪声的影响,Strangfeld对其实验结果进行了2阶Fourier模态的拟合, 图中点线为Strangfeld拟合的结果。

图4为升力系数Cl。在2°和8°迎角时, 升力系数随迎角增长为线性关系。计算所得到的结果振幅略高于实验和理论, 在速度减小时升力系数减小相对滞后。这是由于在马赫数0.0326时流动具有弱可压缩性, 入口速度变化需要一定时间传递到翼型附近。

图5为力矩系数Cm。力矩关于翼尖取矩。与图4计算得到的升力系数情况相似。计算得到的力矩系数的振幅略高于模型实验相近, 但是相位落后。在马赫数0.0326时, 计算得到的非定常升力系数Cl、力矩系数Cm略高于模型和实验, 相位落后, 但误差小于5%, 满足对工程问题的分析的精度要求。下面将采用相同的计算条件, 在马赫数提高为0.1、0.2和0.3时, 分析马赫数提高对升力系数的影响。

马赫数小于0.3时, 流动近似为不可压缩流动。在定常状态下, 马赫数为0.1、0.2、0.3的升力系数与马赫数为0.0326时相同。图6为来流速度正弦振动时升力系数随时间的变化。图中比较了马赫数0.0326、0.1、0.2、0.3时非定常升力系数与Isaacs、Greenberg模型的差别。图6(a)的迎角为2°; 图6(b)的迎角为8°。由于当马赫数升高时, 实验的条件不再适用, 所以不将实验与计算进行比较。随着马赫数的升高, 非定常升力系数的最高点逐渐升高, 相位逐渐落后, 其中对于最大马赫数0.3时与模型的差别最大。

马赫数变化, 同样会引起力矩系数的变化。图7为来流速度正弦振动时升力系数随时间的变化。随着马赫数的升高, 力矩系数与模型的差别逐渐增加。其变化趋势与升力系数相同。在最大马赫数0.3时差别对大。对于相同的马赫数, 差别最明显的位置位于速度最大的点, 即ωt=π/2。而对于速度最小的时候,ωt=3π/2时, 速度模型与计算得到的力矩系数最相近。

由于计算与模型考虑的是非定常情况下的升力系数、力矩系数采用了定常情况下的升力系数、力矩系数作归一化。所以图6与图7中马赫数引起的变化即可能为分子——非定常时升力、力矩系数的变化, 也可能是分母——定常升力、力矩系数的变化。图8为在定常情况下将马赫数为0.1、0.2和0.3的升力与阻力系数与马赫数为0.0326的情况相比。图中方形表示升力系数之比, 三角形表示力矩系数之比。实心图形表示迎角为2°的情形, 空心图形表示迎角为8°的情况。

从图中可以看出, 随着马赫数的增大, 升力系数与力矩系数逐渐增大, 与马赫数为0.0326时的差别也逐渐扩大。但是误差在5%以内, 可以认为近似相等。也就是升力系数与力矩系数在定常情况下随着马赫数的增加基本保持不变。图6和图7中不同马赫数时升力系数与力矩系数的差别主要由于非定常部分引起。

根Isaacs的理论[2], 非定常流动产生的升力与环量和速度相关。

L(ωt)=ρu(ωt)Γ(ωt)+

式中L(ωt)为非定常升力,Γ为绕翼型的环量,γb为受约束涡层(boundvorticitysheet),c为翼型的弦长。Isaacs的模型与Greenberg的模型考虑了速度与环量受到的非定常效应的影响。不可压缩流动假设下密度近似为常数不随压力的变化而改变。在本文计算设置下, 入口速度的变化会转化为压力的变化从而驱动流场中速度发生变化。当考虑流体的可压缩性时, 压力变化会引起密度的变化。相对翼型, 入口来流速度变化时, 也引起了来流的密度发生了变化。下面将分析在入口速度发生变化时, (10)式中密度与来流马赫数之间的关系, 并对Isaacs的模型与Greenberg的模型进行马赫数变化引起密度变化的修正, 使其能够反映来流马赫数对非定常升力系数的影响。

3.2 弱可压缩性修正

假设, 弱可压缩与不可压缩的速度与涡量场相同, 可以采用Isaacs或Greenberg模型描述。弱可压缩性仅改变了式(10)中的密度项。设不可压缩流动假设下模型预测的升力为LI, 弱可压缩假设下得到的升力为L, 根据(10)式得到两者之比R(ωt)为密度之比:

式中ρ为弱可压缩时得到的流动密度, 不可压缩假设下流动的密度为来流密度ρ∞。

将密度ρ分解为两部分, 一部分为来流密度ρ∞, 一部分为流动引起的密度的变化ρ′, 即ρ=ρ∞+ρ′。将压力也分解为两部分, 一部分为来流压力p∞, 一部分为流动引起的压力变化p′, 即p=p∞+p′。在势流假设下流动沿流线满足Bernoulli定理。当马赫数小于0.3时, 定常流动近似为不可压缩。此时的压力脉动量级为:

根据式(12), 流场中的温度假设为常数。对于等温情况, 压力变化与来流压力之比等于密度变化与来流密度的比值:

将等式(15)两边同时减去1得到。

将等温过程压力与密度关系式(16), 代入式(14)得到:

化简得到:

将式(18)等式两边同时加1, 并代入式(11), 得到的压力与模型中的压力的比值为:

由于升力系数为升力与参考面积、来流动压之比。弱可压缩假设与不可压缩假设的参考面积、来流动压相同。所以根据式(19)Isaacs、Greenberg模型预测的升力系数Cl,I与计算得到的升力系数Cl有关系式:

力矩为压力与矩心距离乘积的积分, 不可压缩与弱可压缩情况力矩系数采用的无量纲化系数相同, 因此力矩系数之比等于压力之比。根据式(15)与式(11), 不可压缩与可压缩时的压力之比等于升力之比R(ωt)。

所以力矩之比为:

将迎角为2°和8°时马赫数为0.0326、0.1、0.2和0.3情况下得到的力矩系数除以式(20)R(ωt)得到图10。R(ωt)中的系数K同样取为K=1.8。与升力系数情况相同, 修正使得力矩系数变化的幅值与模型相差在5%以内。

4 结 论

本文采用计算方法, 针对NACA0018翼型, 在来流马赫数分别为0.0326、0.1、0.2和0.3时得到了非定常升力系数随时间的演化。发现:

1) 在马赫数0.0326时计算结果与Isaacs的模型、Greenberg的模型以及Strangfeld等[12]在2014年实验的误差在5%以内。

2) 随着马赫数的升高, 计算与Isaacs的模型、Greenberg的模型的差别逐渐增加, 在马赫数0.3时差别最大, 最大相差50%。

3) 通过考虑入口附近速度变化引起的密度变化, 修正翼型的来流密度, 使得计算结果与Isaacs的模型、Greenberg的模型差别在5%左右。

致谢:感谢西北工业大学王刚教授、刘毅博士在本文工作中给予作者的帮助与建议。

[1]Theodorsen T. General theory of aerodynamic instability and the mechanism of flutter[R]. NACA Rep No. 496, 1935: 413-433.

[2]Isaacs R. Airfoil theory for flows of variable velocity[J]. Journal of the Aeronautical Sciences, 1945, 12(1): 113-117.

[3]Isaacs R. Airfoil theory for rotary wing aircraft[J]. Journal of the Aeronautical Sciences, 1946, 13(4): 218-220.

[4]Greenberg J M. Airfoil in sinusoidal motion in a pulsating stream[R]. NACA TN1326, 1947.

[5]Mateescu D, Neculita S. Low frequency oscillations of thin airfoils in subsonic compressible flows[C]//44th AIAA Aerospace Sciences Meeting and Exhibit, Reno, Nevada, 2006.

[6]William P Walker, Mayuresh J Patil. Unsteady aerodynamics of deformable thin airfoils[J]. Journal of Aircraft, 2014, 51

(6): 1673-1680.

[7]Ramesh K, Gopalarathnam A, Edwards J R, et al. Theoretical analysis of perching and hovering maneuvers[R]. AIAA 2013-3194.

[8]Cho S H, Kim T, Song S J. Freestream Pulsation Effects on the Aeroelastic Response of a Finite Wing[J]. AIAA Journal, 2008, 46(11): 2723-2729.

[9]刘启宽, 李亮, 张志军, 等. 风力机叶片大挠度挥舞振动特性分析[J]. 动力学与控制学报, 2012, 10(2): 171-177.Liu Q, Li L, Zhang Z, et al. Flapwise characteristics of a wind turbine blade with large deflection[J]. Journal of Dynamics and Control, 2012, 10(2): 171-177.

[10]Williams D, Quach V, Kerstens W, et al. Low-Reynolds number wing response to an oscillating freestream with and without feed forward control[R]. AIAA 2009-143.

[11]Jose A I, Baeder J D. Steady and unsteady aerodynamic modeling of trailing edge flaps with overhang and gap using CFD and lower order models[R]. AIAA 2009-1071.

[12]Strangfeld C, Müller-Vahl H, Greenblatt D, et al. Airfoil subjected to high-amplitude free-stream oscillations: theory and experiments [C]//7th AIAA Theoretical Fluid Mechanics Conference, Atlanta, GA, 2014.

[13]Mian H H, Gang Wang, Raza M A. Application and validation of HUNS3D flow solver for aerodynamic drag prediction cases[C]//Proceedings of 2013 10th International Bhurban Conference on Applied Sciences and Technology, Islamabad, Pakistan, 2013: 209-218.

[14]Liou Meng Sing. A sequel to AUSM, Part II: AUSM+-up for all speeds[J]. Journal of Computational Physics, 2006, 214: 137-170.

[15]Spalart P R, Allmaras S R. A one-equation turbulence model for Aerodynamic flows[J]. Le Recherche Aerospatiale, 1994, 1: 5-21.

[16]Katz J, Plotkin A. Low-speed aerodynamics[M]. Cambridge, Cambridge University Press, 2001.

CFD verification and weak compressibility correction of unsteady aerodynamic force models applied to high-amplitude oscillating incoming flows

Guo Li, Lyu Jinan, Feng Feng, Wang Qiang*

(ChinaAcademyofAerospaceAerodynamics,Beijing100074,China)

The two-dimensional airfoil theories of Isaacs and Greenberg for unsteady aerodynamic forces are widely adopted to estimate the aerodynamic performance of wind-turban blades and helicopter blade loads. The models are established under the assumption that the incoming flow is incompressible and without viscosity. However, the viscosity and compressibility are inevitable and the applicability of the model to predict the aerodynamic force in real flows needs to be checked. For the viscous effects, Strangfeld et al. verified the models experimentally using the data of NACA 0018 from wind tunnel at Reynolds number 0.25 million in 2014. The Mach number of the experiment is near 0.0326, which makes the flow almost incompressible. To check the effects of compressibility, a numerical simulation of NACA 0018 is conducted. For verification, the result of Strangfeld et al. at Mach number 0.0326 is repeated using CFD. The simulation further extends to the Mach number 0.1, 0.2 and 0.3 cases to investigate performance of the model at higher Mach numbers. The results show that maximum lift coefficient increases and the phase lags with Mach number increasing. Maximum overshot is at Mach number 0.3, which is over 50% more than the prediction of the models. The results show that the compressibility is non-negligible and the Mach number is needed to be taken into account as a sensitive factor. To extend the application range of Mach number, a correction is added to the models to account the influence of weak compressibility. The correction relates the variation of density of the incoming flow to the Mach number changing, which changes the density near the airfoil. As the aerodynamic force are proportional to the density, the correction would increase the aerodynamic force as Mach number increases. Using the correction, the differences of the models and simulation is within 5%.

Isaacs′s model; Greenberg′s model; unsteady aerodynamic force; low Mach number; compressibility correction

0258-1825(2017)01-0093-08

2015-05-25;

2015-10-26

国家自然科学重大研究计划重点项目(91216202), 国家高技术研究发展计划(863计划) (2015AA01A302)

郭力(1984-), 男, 博士, 工程师. 研究方向: 气动弹性数值模拟. E-mail: gargo@sina.com

王强*(1967-), 男, 博士, 研究员. E-mail: qwang327@163.com

郭力, 吕计男, 冯峰, 等. 大振幅振荡来流条件下非定常气动力模型计算验证与弱可压缩性修正[J]. 空气动力学学报, 2017, 35(1): 93-100.

10.7638/kqdlxxb-2015.0065 Guo L, Lyu J N, Feng F, et al. CFD verification and weak compressibility correction of unsteady aerodynamic force models applied to high-amplitude oscillating incoming flows[J]. Acta Aerodynamica Sinica, 2017, 35(1): 93-100.

V211.3

A doi: 10.7638/kqdlxxb-2015.0065