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

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

?

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

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

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

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

關(guān)鍵詞:俯仰動導(dǎo)數(shù);時間步長;內(nèi)迭代次數(shù);強迫運動振幅;減縮頻率

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

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

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

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

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

(1)

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

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

(2)

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

(3)

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

(4)

當(dāng)非定常運動時間足夠長時,非定常氣動力矩會表現(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為一個運動周期內(nèi)的迭代次數(shù)。

(8)

1.2數(shù)值方法

笛卡爾坐標(biāo)下,三維N-S方程可以寫為

(9)

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

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

2建模與定常計算結(jié)果討論

2.1幾何模型及網(wǎng)格劃分

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

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

圖1　計算模型

圖2　流場網(wǎng)格

2.2定常計算結(jié)果

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

4標(biāo)模俯仰動導(dǎo)數(shù)計算結(jié)果與實驗結(jié)果比較

采用上述方法,取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條件下標(biāo)模的動導(dǎo)數(shù)進行計算,CFD計算結(jié)果、實驗結(jié)果[6]及文獻[9]計算結(jié)果如圖9所示?？芍?計算結(jié)果精度較高,Ma≥1.5時,與實驗結(jié)果相比最大相對誤差不超過5%,與文獻計算結(jié)果基本重合。0.9≤Ma<1.5時,實驗結(jié)果存在較大的散布,本文計算結(jié)果與實驗結(jié)果的平均值較為接近。Ma<0.9時,本文計算結(jié)果與實驗結(jié)果吻合較好,優(yōu)于文獻[9]計算結(jié)果。

圖9　俯仰阻尼動導(dǎo)數(shù)結(jié)果比較

5結(jié)束語

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

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

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

參考文獻

[1]袁先旭,張涵信.基于CFD方法的俯仰靜/動導(dǎo)數(shù)數(shù)值計算[J].空氣動力學(xué)學(xué)報,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的動導(dǎo)數(shù)計算與減縮頻率影響分析[J].飛行力學(xué),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]劉緒.高超聲速內(nèi)外流一體化飛行器動態(tài)特性研究[D].長沙:國防科學(xué)技術(shù)大學(xué),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]盧學(xué)成.超音速、高超音速飛行器動導(dǎo)數(shù)高效的計算方法[J].航空計算技術(shù),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- ),男,工程師,研究方向為氣動設(shè)計與分析。E-mail:wubin_no8@qq.com。

中圖分類號:V211.3

文獻標(biāo)識碼:A

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