伍　彬,陸　韻,周志超,傅建明

(上海機電工程研究所,上海 201100)

?

基于非定常CFD的俯仰動導數(shù)計算方法

伍彬,陸韻,周志超,傅建明

(上海機電工程研究所,上海 201100)

摘要:采用基于N-S方程求解的非定常CFD技術,對俯仰動導數(shù)數(shù)值計算方法進行了研究。研究了時間步長、內迭代次數(shù)、強迫運動幅值、減縮頻率對俯仰動導數(shù)計算結果的影響。計算了動導數(shù)標模在不同馬赫數(shù)下的俯仰動導數(shù),并與實驗結果進行了比較。結果表明,時間步長與減縮頻率是影響俯仰動導數(shù)計算結果的2個重要因素。給出了俯仰動導數(shù)數(shù)值計算時間步長與減縮頻率選取的一般原則?；诜嵌ǔFD技術的俯仰動導數(shù)數(shù)值計算方法精度較高,與實驗結果吻合良好,具有較高的工程應用價值。

關鍵詞:俯仰動導數(shù);時間步長;內迭代次數(shù);強迫運動振幅;減縮頻率

動導數(shù)是表征飛行器動態(tài)特性的重要參數(shù),對飛行器的穩(wěn)定性和操縱性有著十分重要的影響。目前,獲取動導數(shù)的主要方法包括理論估算、風洞實驗、飛行實驗及數(shù)值計算。其中理論估算最為快速方便,但只適用于簡單、常規(guī)的飛行器,且精度較差,難以滿足復雜的現(xiàn)代飛行器的需求。風洞實驗可以通過有限自由振動法、強迫運動法和模型自由飛等手段獲得飛行器動導數(shù),但面臨耗費大、周期長的問題。而飛行實驗難度大、周期長、風險高,且不能在飛行器初步設計階段獲得數(shù)據,無法對設計起指導作用。隨著計算流體力學(CFD,Computational Fluid Dynamics)技術的發(fā)展及非定常數(shù)值計算精度的提高,采用CFD數(shù)值計算的方法獲取飛行器動導數(shù)是目前較理想、較經濟的方法之一。

袁先旭、張涵信[1]等研究了基于CFD的俯仰靜、動導數(shù)計算,提出了2種辨識方法。孫濤、高正紅[2]等基于Euler方程研究了動導數(shù)的數(shù)值計算方法,以及Ma為1.58的情況下不同減縮頻率對動導數(shù)計算結果的影響。劉緒[3]、盧學成[4]分別研究了超聲速情況下的飛行器動導數(shù)數(shù)值計算。Despirito[5]和Bhagwandin[6]分別采用錐動法對動導數(shù)進行了計算,得到了較好的結果。本文基于N-S方程對飛行器的動導數(shù)計算方法開展了研究,以動導數(shù)標模[7-8](Basic Finner Missile)為研究對象,采用小振幅強迫運動的方法獲取飛行器的俯仰動導數(shù),研究了時間步長、內迭代次數(shù)、強迫運動振幅以及減縮頻率對動導數(shù)計算結果的影響,給出了一般性的結論。計算了動導數(shù)標模在不同馬赫數(shù)下的俯仰動導數(shù),并與實驗結果進行了比較。

1動導數(shù)計算理論與數(shù)值方法

1.1動導數(shù)計算理論

以俯仰阻尼動導數(shù)計算為例,對動導數(shù)的計算方法進行了說明。假設飛行器在俯仰平面內做繞質心的單自由度小振幅強迫運動。強迫運動的具體形式為

(1)

式中:θ為俯仰角,θm為振幅,α為攻角,ω為強迫運動頻率,q為俯仰角速度。

飛行器所受的非定常俯仰力矩系數(shù)Cm的表達式可以寫成[3]:

(2)

將式(2)進行Taylor級數(shù)展開并略去高階項,有:

(3)

將式(1)代入式(3),合并同類項,有:

(4)

當非定常運動時間足夠長時,非定常氣動力矩會表現(xiàn)出穩(wěn)定的周期性,且其周期與強迫運動周期相同,因此非定常俯仰力矩系數(shù)也可以寫為

Cm=Cm0+Cmssinωt+Cmccosωt

(5)

式中:Cms為正弦項系數(shù),Cmc為余弦項系數(shù)。

對比式(4)和式(5),可知:

(6)

非定常俯仰力矩系數(shù)隨時間的變化規(guī)律可以由CFD計算得到,那么Cms和Cmc可以由下式給出:

(7)

式中:N為一個運動周期內的迭代次數(shù)。

(8)

1.2數(shù)值方法

笛卡爾坐標下,三維N-S方程可以寫為

(9)

式中:Q為守恒量;F1,F2,F3為無粘通量;G1,G2,G3為粘性通量。

本文計算時,空間離散采用二階Roe格式,時間推進采用隱式雙時間步推進方法,在時間上也保證二階精度。湍流模型采用兩方程Realizablek-ε模型,該湍流模型能夠較好地處理復雜剪切流、漩渦流動、邊界層分離、大攻角失速和轉捩等狀況。強迫運動采用剛性動網格技術實現(xiàn)。

2建模與定常計算結果討論

2.1幾何模型及網格劃分

計算模型采用動導數(shù)標模,彈徑為d,長細比為10。尾翼厚度為0.08d,尾翼前緣倒圓半徑為0.004d。具體尺寸見圖1。

劃分流場網格,第一層網格高度約為0.005 mm,網格總數(shù)為280萬,物面及對稱面網格如圖2所示。

圖1　計算模型

圖2　流場網格

2.2定常計算結果

在非定常計算之前,先進行定常計算,并與實驗結果進行比較,以確認定常數(shù)值計算的精度,并將定常流場作為非定常計算的初始流場。圖3給出了Ma=0.9和Ma=2.5兩個狀態(tài)下的對稱面流場馬赫數(shù)分布云圖和等值線圖,可以看出,頭部激波和底部低速區(qū)等典型流動特征都得到了很好的捕捉。

圖3　馬赫數(shù)等值線及物面壓力分布

圖4給出了Ma=0.6～3.0范圍內飛行器的軸向力系數(shù)CA、法向力系數(shù)斜率CNα和俯仰力矩系數(shù)斜率Cmα。從圖4可以看出,定常計算的結果與實驗結果吻合很好,具有足夠的精度,能滿足非定常計算的要求。

圖4　定常氣動系數(shù)計算值與實驗結果比較

3動導數(shù)計算影響因素分析

影響動導數(shù)計算結果的因素有網格類型、網格數(shù)量、湍流模型、非定常時間步長、內迭代次數(shù)、強迫運動幅值以及減縮頻率等,本文只考慮非定常時間步長、內迭代次數(shù)、強迫運動幅值以及減縮頻率的影響,對動導數(shù)的影響因素進行分析。

3.1時間步長與內迭代次數(shù)的影響

定義Δt為非定常計算時間步長,NI為內迭代次數(shù)。圖5給出了不同時間步長的俯仰力矩系數(shù)隨攻角的變化結果,圖6給出了不同時間步長和內迭代次數(shù)的俯仰動導數(shù)計算結果?？芍?當強迫運動幅值和減縮頻域一定時,不同時間步長對應的攻角-俯仰力矩系數(shù)遲滯曲線基本相同,但俯仰動導數(shù)計算結果區(qū)別較大。隨著時間步長減小,計算結果相對差別也逐漸減小。Ma=0.9時,相對差別從3%下降到1%;Ma=2.5時,相對差別從2%下降到0.7%。內迭代次數(shù)對結果影響較小,當內迭代次數(shù)NI≥10時,俯仰阻尼動導數(shù)結果趨于穩(wěn)定,不隨內迭代次數(shù)的增加而變化。

圖5　不同時間步長情況下俯仰力矩系數(shù)隨攻角的變化曲線

圖6　不同時間步長和內迭代次數(shù)情況下的俯仰阻尼動導數(shù)

3.2強迫運動幅值與減縮頻率的影響

圖7給出了不同強迫運動幅值的俯仰力矩系數(shù)隨攻角的變化結果,圖8給出了不同強迫運動幅值和減縮頻率的俯仰動導數(shù)計算結果。

圖7　不同強迫運動幅值情況下俯仰力矩系數(shù)隨攻角的變化曲線

圖8　不同運動幅值和減縮頻率情況下的俯仰阻尼動導數(shù)

由圖可知,強迫運動幅值越大,俯仰力矩系數(shù)越大。運動幅值對俯仰動導數(shù)計算結果影響很小,當時間步長和內迭代次數(shù)一定時,不同運動幅值對應的俯仰動導數(shù)數(shù)值計算結果很接近。減縮頻率對動導數(shù)計算結果有著顯著影響,當減縮頻率超出某個范圍時,動導數(shù)計算結果逐漸偏離實驗值,計算精度迅速下降。其原因是減縮頻率k很大時,式(3)略去與k2有關的高階項會導致誤差增大;減縮頻率k很小時,非定常遲滯效應不能準確體現(xiàn),也會導致誤差增大。在實際計算時,應選擇與實驗減縮頻率值相近的減縮頻率進行動導數(shù)計算,才能保證計算結果與實驗結果吻合較好。

4標模俯仰動導數(shù)計算結果與實驗結果比較

采用上述方法,取N=200,NI=20,θm=0.25°,取k=0.01,對Ma=0.6,0.8,0.9,1.05,1.5,2.0,2.5,3.0條件下標模的動導數(shù)進行計算,CFD計算結果、實驗結果[6]及文獻[9]計算結果如圖9所示?？芍?計算結果精度較高,Ma≥1.5時,與實驗結果相比最大相對誤差不超過5%,與文獻計算結果基本重合。0.9≤Ma<1.5時,實驗結果存在較大的散布,本文計算結果與實驗結果的平均值較為接近。Ma<0.9時,本文計算結果與實驗結果吻合較好,優(yōu)于文獻[9]計算結果。

圖9　俯仰阻尼動導數(shù)結果比較

5結束語

本文結合CFD技術對動導數(shù)數(shù)值計算方法進行了研究,討論了時間步長、內迭代次數(shù)、強迫運動幅值、減縮頻率對俯仰動導數(shù)計算結果的影響,計算了標模的俯仰動導數(shù),與實驗結果進行了比較。得到以下結論:

①時間步長與減縮頻率是影響俯仰動導數(shù)計算結果的2個關鍵因素,應結合實驗情況選擇適當?shù)臅r間步長和減縮頻率進行計算,兼顧計算精度與計算效率。

②通過對標模的俯仰動導數(shù)進行計算,并與實驗結果對比,可知本文采用的動導數(shù)數(shù)值計算方法能夠得到較為精確的計算結果,能夠滿足型號研制的需求。

參考文獻

[1]袁先旭,張涵信.基于CFD方法的俯仰靜/動導數(shù)數(shù)值計算[J].空氣動力學學報,2005,23(4):458-463.

YUAN Xian-xu,ZHANG Han-xin.The pitching satic/dynamic derivatives computation based on CFD method[J].Acta Aerodynamics Sinica,2005,23(4):458-463.(in Chinese)

[2]孫濤,高正紅.基于CFD的動導數(shù)計算與減縮頻率影響分析[J].飛行力學,2011,29(4):15-18.

SUN Tao,GAO Zheng-hong.Identify of aircraft dynamic derivatives based on CFD technology and analysis of reduce frequency[J].Flight Dynamics,2011,29(4):15-18.(in Chinese)

[3]劉緒.高超聲速內外流一體化飛行器動態(tài)特性研究[D].長沙:國防科學技術大學,2011:14-17.

LIU Xu.Investigation of dynamic characteristics of hypersonic airframe/propulsion integrative vehicle[D].Changsha:National University of Defense Technology,2011:14-17.(in Chinese)

[4]盧學成.超音速、高超音速飛行器動導數(shù)高效的計算方法[J].航空計算技術,2008,38(3):28-31.

LU Xue-cheng.A high efficient method for computing dynamic derivatives of supersonic/hypersonic aircraft[J].Aeronautical Computing Technique,2008,38(3):28-31.(in Chinese)

[5]DESPIRITO J,SILTON S I,WEINACHT P.Navier-Stokes predictions of dynamic stability derivatives evaluation of steady-state methods[J].Journal of Spacecraft and Rockets,2009,46(6):1-13.

[6]VISUAL A B,JUBARAJ S.Numerical prediction of pitch damping stability derivatives for finned projectiles,AIAA 2011-0328[R].2011.

[7]Arnold Engineering Development Center.Experimental roll-damping,mangus,and static-stability characterictics of two slender missile configurations at high angles of attack(0 to 90 Deg)and Mach numbers 0.2 to 2.5,AD-A027[R].1976.

[8]DUPUIS A D,HATHAWAY W.Aeroballistic range and wind tunnel tests of the basic finner reference projectile from subsonic to high supersonic velocities,TM 2002-136[R].2002.

[9]JUBARAJ S.Numerical computations of dynamic derivatives of a finned projectile using a time-accurate CFD method,AIAA 2007-6581[R].2007

Calculation Method of Pitching Dynamic Derivatives Based on Unsteady CFD Technology

WU Bin,LU Yun,ZHOU Zhi-chao,FU Jian-ming

(Shanghai Electro-mechanical Engineering Institute,Shanghai 201100,China)

Abstract:Based on solving N-S equations,the numerical prediction method of pitching dynamic derivatives was studied.The effects of transient time-step,inner iterations,forced sinusoidal motion amplitude and reduced frequency on the calculation result of pitching dynamic derivatives were analyzed.Prediction of dynamic derivatives was conducted for the standard model of dynamic derivatives.The result shows that the transient time-step and reduced frequency are the most important factor affecting dynamic derivatives prediction.The general principle of choosing transient time-step and reduced frequency was presented.The calculation method of pitching dynamic derivatives based on unsteady CFD has high accuracy,and the result accords with the experiment data.The method has practical value.

Key words:pitching dynamic derivative;time step;inner iteration;forced sinusoidal motion amplitude;reduced frequency

收稿日期:2016-03-25

作者簡介:伍彬(1987- ),男,工程師,研究方向為氣動設計與分析。E-mail:wubin_no8@qq.com。

中圖分類號:V211.3

文獻標識碼:A

文章編號:1004-499X(2016)02-0042-05