大振幅振蕩來流條件下非定常氣動力模型計算驗證與弱可壓縮性修正

郭 力, 呂計男, 馮 峰, 王 強

(中國航天空氣動力技術(shù)研究院, 北京 100074)

大振幅振蕩來流條件下非定常氣動力模型計算驗證與弱可壓縮性修正

郭 力, 呂計男, 馮 峰, 王 強*

(中國航天空氣動力技術(shù)研究院, 北京 100074)

針對大振幅振蕩來流條件下薄翼受到的非定常氣動力，Isaccs和Greenberg分別發(fā)展了非定常氣動力模型。這兩種模型可以用于直升飛機槳葉與風(fēng)力發(fā)電葉片的氣動力分析，模型在不可壓縮無黏來流條件下建立，但實際流動不可避免粘性和弱可壓縮性的影響，需要檢驗兩種模型的適用性。針對粘性效應(yīng)的影響，2014年Strangfeld對于NACA0018翼型，通過風(fēng)洞實驗驗證了在Reynolds數(shù)25萬時，Isaccs和Greenberg的模型仍適用，實驗的Mach數(shù)為0.0326，流動近似不可壓縮流動。針對可壓縮性的影響，通過數(shù)值模擬方法進行了研究。首先重復(fù)了實驗在Mach數(shù)為0.0326時的結(jié)果，并進一步考察了當(dāng)Mach數(shù)提高為0.1、0.2和0.3時非定常氣動力的變化。結(jié)果表明隨著Mach數(shù)的提高，升力系數(shù)的最高點逐漸高于模型，并且相位逐漸落后，在Mach數(shù)為0.3時差別最明顯，非定常升力系數(shù)最高點計算與模型相差50%。此即表明弱可壓縮性對模型的預(yù)測結(jié)果影響不可忽略。為了擴展模型在Mach數(shù)變化時的適用范圍，對模型進行了弱可壓縮性修正。通過考慮速度變化引起均勻來流中密度的變化，修正了翼型附近流體密度，使其跟隨來流Mach數(shù)變化。采用此方法，將計算與模型的幅值差別減小到5%左右。

Isaacs模型；Greenberg模型；非定常氣動力；低馬赫數(shù)；弱可壓縮性修正

0 引 言

旋轉(zhuǎn)的葉片, 如直升飛機的旋翼, 風(fēng)力機的葉片等, 即使在定常運動時, 由于旋轉(zhuǎn)的作用, 其截面受到的來流大小與迎角也會非定常變化。針對翼型的平動與非定常迎角變化問題, Theodorsen在1935年提出了非定常氣動力的理論解[1]。Isaacs根據(jù)Theodorsen理論加入了來流速度變化的影響[2]與迎角變化的影響[3]。Greenberg考慮了速度脈動與翼型運動的復(fù)合作用[4]。Mateescu在2006年發(fā)展了在低速可壓縮情況下非定常氣動力模型[5], 模型采用在翼型前緣與脊(ridge)上的速度奇點的解, 考察了剛性與柔性翼型的平動與轉(zhuǎn)動非定常性。Walker等考慮了翼型在柔性變形情況下的氣動力[6]。Ramesh發(fā)展了模型用于考察前緣渦對于昆蟲飛行的影響[7]。Cho等采用3維非定常渦格(vortex lattice)方法, 發(fā)展了非定常氣動力模型, 用于考慮速度脈動引起的非定常氣動力對氣動彈性穩(wěn)定性的影響[8]。劉啟寬等[9]采用Greenberg公式分析了風(fēng)力機葉片的揮舞特性。Williams等采用Greenberg模型[4]設(shè)計控制率, 抑制來流振蕩引起的受力變化[10]。Jose等采用Theodorsen的理論分析了二維多段翼型間隙對非定常氣動力的影響[11]。Strangfeld[12]使用NACA0018翼型通過實驗對Iscaccs和Greenberg提出的模型進行了考察。實驗中翼型固定, 來流速度做正弦變化。來流的平均速度為11m/s, 最大速度為16.5m/s, 最小速度為5.5m/s。在海平面條件下, 平均速度對應(yīng)馬赫數(shù)近似為0.0326, 流動可以近似為不可壓縮流動?；谙议L的Reynolds數(shù)為2.5×105。實驗考察了迎角為2°與8°時升力系數(shù)隨速度的變化。在迎角為2°時實驗結(jié)果振動較大。但通過前兩階Fourier模態(tài)的擬合結(jié)果與理論吻合較好。在迎角為8°時實驗較光滑, 但低估了最大升力系數(shù)。這是由于迎角8°時實驗中渦在翼型表面的分布與理論分布不同引起。

Isaacs和Greenberg的氣動力公式采用勢流理論推導(dǎo)。其試用范圍為: 二維不可壓縮勢流流動, 迎角α較小, 滿足近似關(guān)系α≈sin(α), 來流不會反向。本文通過計算, 考察弱可壓縮性對Isaacs和Greenberg模型精度的影響。文中首先重復(fù)了Strangfeld等的實驗條件下升力系數(shù)隨速度的變化。然后對更高馬赫數(shù)0.1, 0.2和0.3時升力系數(shù)的變化進行考察。

1 Greenberg、Isaacs非定常升力系數(shù)氣動力模型

來流速度圍繞平均速度做正弦振動時會引起翼型受到的力和力矩也隨著速度變化, 做正弦變化。設(shè)來流速度u做正弦變化的表達式為:

Cl、Cm和lm的表達式為無窮級數(shù),Strangfeld等[12]將Cl中的無窮級數(shù)在第8項之后(m>8)截斷;lm的級數(shù)表達式在第15項(n>15)后截斷。其誤差小于10-6。

根據(jù)Greenberg模型, 非定常升力系數(shù)Cl(t)與定常升力系數(shù)Cl,st之比為:

式中a是取矩的點距離前緣的距離。在不可壓縮無粘勢流假設(shè)下, 對于翼型四分之一弦長的地方取矩, 無論迎角多少, 得到的力矩系數(shù)近似為0[16]。文中所有力矩均關(guān)于翼型的前緣取矩。

圖1為Iscaacs和Greenberg模型的比較。其中實線為Isaacs模型結(jié)果, 虛線為Greenberg模型結(jié)果。無論對于升力系數(shù), 還是力矩系數(shù),Greenberg公式在最低點與最高點處低估了升力。但是在初始位置與接近平衡點時兩個模型結(jié)果相近。

2 翼型繞流計算設(shè)置

Isaacs模型與Greenberg模型均對于不可壓縮流動提出。為了與模型的適用范圍接近, 本文計算采用的馬赫數(shù)較低。分別考察了馬赫數(shù)為0.0326、0.1、0.2和0.3時, 迎角為2°和8°的情況。其中馬赫數(shù)為0.0326近似為Strangfeld實驗[12]的來流條件。計算采用HUNS3D程序[13]。通量部分采用AUSM+-up格式[14], AUSM+-up格式在馬赫數(shù)較低時可以保證格式的收斂性和精度[14]。湍流模型采用SA模型[15], Reynolds數(shù)為Re=2.5×105。時間推進采用雙時間步方法, 其中隱式迭代采用LU-SGS方法。計算采用2維網(wǎng)格, 計算區(qū)域為矩形, 如圖2所示。

在流向與法向方向上為40倍的翼型弦長, 出口的上下兩角導(dǎo)圓。翼型位于計算區(qū)域中央。迎角通過翼型旋轉(zhuǎn)設(shè)定。來流方向始終與入口垂直, 上下兩個邊界采用滑移無穿透邊界條件, 使得入口下游的流場能夠根據(jù)來流速度的變化進行調(diào)整, 避免人為指定帶來的誤差。入口邊界指定速度, 其隨時間的變化為式(1)。式(1)中σ=0.5與Strangfeld等的實驗[12]所用振幅相同,ω=4.72×10-3根據(jù)Strangfeld的實驗[12]中約化頻率k=0.0074得到。

為了驗證網(wǎng)格與計算條件的適用性, 采用程序?qū)砹黢R赫數(shù)為0.0326、迎角為8°的定常情況進行計算, 并與實驗進行對比。圖3為翼型表面升力系數(shù)。其中實線為計算所得結(jié)果, 點為實驗結(jié)果。在翼型上表面與下表面上計算與實驗吻合較好。在前緣處, 下表面得到的升力系數(shù)與實驗接近。在前緣上表面和尾緣處壓力略高于實驗所得值, 但差別在誤差范圍內(nèi)。

3 非定常氣動力特性分析

3.1 模型與計算比較

為了驗證程序在非定常時的計算結(jié)果, 將計算得到的升力系數(shù)Cl、力矩系數(shù)Cm在馬赫數(shù)0.0326時與實驗和理論進行了比較。為了方便比較, 升力系數(shù)采用定常狀態(tài)做了歸一化。非定常計算從定常狀態(tài)啟動, 為了消除啟動時初始效應(yīng)對計算的影響, 結(jié)果為速度變化第二個周期時的數(shù)據(jù)。圖4、圖5中實線為Isaacs的理論結(jié)果, 虛線為Greenberg的理論結(jié)果，點劃線為計算結(jié)果, 實心圓點為Strangfeld等的實驗結(jié)果[12]。實驗中壓力傳感器受到噪聲的影響, 出現(xiàn)脈動[12]。為了減小實驗時壓力傳感器噪聲的影響,Strangfeld對其實驗結(jié)果進行了2階Fourier模態(tài)的擬合, 圖中點線為Strangfeld擬合的結(jié)果。

圖4為升力系數(shù)Cl。在2°和8°迎角時, 升力系數(shù)隨迎角增長為線性關(guān)系。計算所得到的結(jié)果振幅略高于實驗和理論, 在速度減小時升力系數(shù)減小相對滯后。這是由于在馬赫數(shù)0.0326時流動具有弱可壓縮性, 入口速度變化需要一定時間傳遞到翼型附近。

圖5為力矩系數(shù)Cm。力矩關(guān)于翼尖取矩。與圖4計算得到的升力系數(shù)情況相似。計算得到的力矩系數(shù)的振幅略高于模型實驗相近, 但是相位落后。在馬赫數(shù)0.0326時, 計算得到的非定常升力系數(shù)Cl、力矩系數(shù)Cm略高于模型和實驗, 相位落后, 但誤差小于5%, 滿足對工程問題的分析的精度要求。下面將采用相同的計算條件, 在馬赫數(shù)提高為0.1、0.2和0.3時, 分析馬赫數(shù)提高對升力系數(shù)的影響。

馬赫數(shù)小于0.3時, 流動近似為不可壓縮流動。在定常狀態(tài)下, 馬赫數(shù)為0.1、0.2、0.3的升力系數(shù)與馬赫數(shù)為0.0326時相同。圖6為來流速度正弦振動時升力系數(shù)隨時間的變化。圖中比較了馬赫數(shù)0.0326、0.1、0.2、0.3時非定常升力系數(shù)與Isaacs、Greenberg模型的差別。圖6(a)的迎角為2°; 圖6(b)的迎角為8°。由于當(dāng)馬赫數(shù)升高時, 實驗的條件不再適用, 所以不將實驗與計算進行比較。隨著馬赫數(shù)的升高, 非定常升力系數(shù)的最高點逐漸升高, 相位逐漸落后, 其中對于最大馬赫數(shù)0.3時與模型的差別最大。

馬赫數(shù)變化, 同樣會引起力矩系數(shù)的變化。圖7為來流速度正弦振動時升力系數(shù)隨時間的變化。隨著馬赫數(shù)的升高, 力矩系數(shù)與模型的差別逐漸增加。其變化趨勢與升力系數(shù)相同。在最大馬赫數(shù)0.3時差別對大。對于相同的馬赫數(shù), 差別最明顯的位置位于速度最大的點, 即ωt=π/2。而對于速度最小的時候,ωt=3π/2時, 速度模型與計算得到的力矩系數(shù)最相近。

由于計算與模型考慮的是非定常情況下的升力系數(shù)、力矩系數(shù)采用了定常情況下的升力系數(shù)、力矩系數(shù)作歸一化。所以圖6與圖7中馬赫數(shù)引起的變化即可能為分子——非定常時升力、力矩系數(shù)的變化, 也可能是分母——定常升力、力矩系數(shù)的變化。圖8為在定常情況下將馬赫數(shù)為0.1、0.2和0.3的升力與阻力系數(shù)與馬赫數(shù)為0.0326的情況相比。圖中方形表示升力系數(shù)之比, 三角形表示力矩系數(shù)之比。實心圖形表示迎角為2°的情形, 空心圖形表示迎角為8°的情況。

從圖中可以看出, 隨著馬赫數(shù)的增大, 升力系數(shù)與力矩系數(shù)逐漸增大, 與馬赫數(shù)為0.0326時的差別也逐漸擴大。但是誤差在5%以內(nèi), 可以認為近似相等。也就是升力系數(shù)與力矩系數(shù)在定常情況下隨著馬赫數(shù)的增加基本保持不變。圖6和圖7中不同馬赫數(shù)時升力系數(shù)與力矩系數(shù)的差別主要由于非定常部分引起。

根Isaacs的理論[2], 非定常流動產(chǎn)生的升力與環(huán)量和速度相關(guān)。

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

式中L(ωt)為非定常升力,Γ為繞翼型的環(huán)量,γb為受約束渦層(boundvorticitysheet),c為翼型的弦長。Isaacs的模型與Greenberg的模型考慮了速度與環(huán)量受到的非定常效應(yīng)的影響。不可壓縮流動假設(shè)下密度近似為常數(shù)不隨壓力的變化而改變。在本文計算設(shè)置下, 入口速度的變化會轉(zhuǎn)化為壓力的變化從而驅(qū)動流場中速度發(fā)生變化。當(dāng)考慮流體的可壓縮性時, 壓力變化會引起密度的變化。相對翼型, 入口來流速度變化時, 也引起了來流的密度發(fā)生了變化。下面將分析在入口速度發(fā)生變化時, (10)式中密度與來流馬赫數(shù)之間的關(guān)系, 并對Isaacs的模型與Greenberg的模型進行馬赫數(shù)變化引起密度變化的修正, 使其能夠反映來流馬赫數(shù)對非定常升力系數(shù)的影響。

3.2 弱可壓縮性修正

假設(shè), 弱可壓縮與不可壓縮的速度與渦量場相同, 可以采用Isaacs或Greenberg模型描述。弱可壓縮性僅改變了式(10)中的密度項。設(shè)不可壓縮流動假設(shè)下模型預(yù)測的升力為LI, 弱可壓縮假設(shè)下得到的升力為L, 根據(jù)(10)式得到兩者之比R(ωt)為密度之比:

式中ρ為弱可壓縮時得到的流動密度, 不可壓縮假設(shè)下流動的密度為來流密度ρ∞。

將密度ρ分解為兩部分, 一部分為來流密度ρ∞, 一部分為流動引起的密度的變化ρ′, 即ρ=ρ∞+ρ′。將壓力也分解為兩部分, 一部分為來流壓力p∞, 一部分為流動引起的壓力變化p′, 即p=p∞+p′。在勢流假設(shè)下流動沿流線滿足Bernoulli定理。當(dāng)馬赫數(shù)小于0.3時, 定常流動近似為不可壓縮。此時的壓力脈動量級為:

根據(jù)式(12), 流場中的溫度假設(shè)為常數(shù)。對于等溫情況, 壓力變化與來流壓力之比等于密度變化與來流密度的比值:

將等式(15)兩邊同時減去1得到。

將等溫過程壓力與密度關(guān)系式(16), 代入式(14)得到:

化簡得到:

將式(18)等式兩邊同時加1, 并代入式(11), 得到的壓力與模型中的壓力的比值為:

由于升力系數(shù)為升力與參考面積、來流動壓之比。弱可壓縮假設(shè)與不可壓縮假設(shè)的參考面積、來流動壓相同。所以根據(jù)式(19)Isaacs、Greenberg模型預(yù)測的升力系數(shù)Cl,I與計算得到的升力系數(shù)Cl有關(guān)系式:

力矩為壓力與矩心距離乘積的積分, 不可壓縮與弱可壓縮情況力矩系數(shù)采用的無量綱化系數(shù)相同, 因此力矩系數(shù)之比等于壓力之比。根據(jù)式(15)與式(11), 不可壓縮與可壓縮時的壓力之比等于升力之比R(ωt)。

所以力矩之比為:

將迎角為2°和8°時馬赫數(shù)為0.0326、0.1、0.2和0.3情況下得到的力矩系數(shù)除以式(20)R(ωt)得到圖10。R(ωt)中的系數(shù)K同樣取為K=1.8。與升力系數(shù)情況相同, 修正使得力矩系數(shù)變化的幅值與模型相差在5%以內(nèi)。

4 結(jié) 論

本文采用計算方法, 針對NACA0018翼型, 在來流馬赫數(shù)分別為0.0326、0.1、0.2和0.3時得到了非定常升力系數(shù)隨時間的演化。發(fā)現(xiàn):

1) 在馬赫數(shù)0.0326時計算結(jié)果與Isaacs的模型、Greenberg的模型以及Strangfeld等[12]在2014年實驗的誤差在5%以內(nèi)。

2) 隨著馬赫數(shù)的升高, 計算與Isaacs的模型、Greenberg的模型的差別逐漸增加, 在馬赫數(shù)0.3時差別最大, 最大相差50%。

3) 通過考慮入口附近速度變化引起的密度變化, 修正翼型的來流密度, 使得計算結(jié)果與Isaacs的模型、Greenberg的模型差別在5%左右。

致謝:感謝西北工業(yè)大學(xué)王剛教授、劉毅博士在本文工作中給予作者的幫助與建議。

[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]劉啟寬, 李亮, 張志軍, 等. 風(fēng)力機葉片大撓度揮舞振動特性分析[J]. 動力學(xué)與控制學(xué)報, 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

國家自然科學(xué)重大研究計劃重點項目(91216202), 國家高技術(shù)研究發(fā)展計劃(863計劃) (2015AA01A302)

郭力(1984-), 男, 博士, 工程師. 研究方向: 氣動彈性數(shù)值模擬. E-mail: gargo@sina.com

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

郭力, 呂計男, 馮峰, 等. 大振幅振蕩來流條件下非定常氣動力模型計算驗證與弱可壓縮性修正[J]. 空氣動力學(xué)學(xué)報, 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