高攀 GAO Pan

（中國特種飛行器研究所，荊門 448035）

0 引言

動(dòng)導(dǎo)數(shù)是飛行器動(dòng)態(tài)特性分析的必要數(shù)據(jù)，是飛行品質(zhì)分析的重要依據(jù)。目前對(duì)動(dòng)導(dǎo)數(shù)獲取方式主要包括飛行試驗(yàn)、風(fēng)洞試驗(yàn)、工程估算、數(shù)值計(jì)算等方法[1]。飛艇與飛機(jī)氣動(dòng)外形和飛行原理都有較大差別，針對(duì)飛機(jī)的動(dòng)導(dǎo)數(shù)工程估算方法并不適用于飛艇；飛行試驗(yàn)和風(fēng)洞試驗(yàn)代價(jià)高昂、試驗(yàn)周期長(zhǎng)；CFD方法與飛行試驗(yàn)、風(fēng)洞試驗(yàn)相比，具有經(jīng)濟(jì)性高、無支架干擾等優(yōu)點(diǎn)，利用CFD求解飛艇動(dòng)導(dǎo)數(shù)已經(jīng)成為一個(gè)重要研究方向。國內(nèi)對(duì)飛行器動(dòng)導(dǎo)數(shù)計(jì)算開展了大量的研究，陶洋、范召林[2]等采用求解ALE形式RANS方程計(jì)算了帶控制舵導(dǎo)彈俯仰和滾轉(zhuǎn)動(dòng)導(dǎo)數(shù)；伍彬、陸韻[3]等采用基于N-S方程的數(shù)值模擬技術(shù)對(duì)俯仰動(dòng)導(dǎo)數(shù)計(jì)算方法進(jìn)行了研究，研究了時(shí)間步長(zhǎng)、內(nèi)迭代次數(shù)、強(qiáng)迫運(yùn)動(dòng)幅值、減縮頻率對(duì)俯仰動(dòng)導(dǎo)數(shù)計(jì)算結(jié)果的影響；孫濤、高正紅[4]等重點(diǎn)分析了減縮頻率對(duì)動(dòng)導(dǎo)數(shù)計(jì)算的影響，提出了利用CFD進(jìn)行動(dòng)導(dǎo)數(shù)計(jì)算時(shí)減縮頻率的選擇原則。國外也對(duì)飛行器動(dòng)導(dǎo)數(shù)數(shù)值計(jì)算方法開展了很多研究[5-6]。本文首先采用俯仰振蕩法對(duì)國際標(biāo)模Finner導(dǎo)彈進(jìn)行動(dòng)導(dǎo)數(shù)計(jì)算，并將計(jì)算結(jié)果與風(fēng)洞試驗(yàn)結(jié)果進(jìn)行對(duì)比，驗(yàn)證動(dòng)導(dǎo)數(shù)CFD求解的準(zhǔn)確性；然后針對(duì)某飛艇開展CFD計(jì)算，采用俯仰振蕩法和旋轉(zhuǎn)流場(chǎng)法求解縱向單獨(dú)動(dòng)導(dǎo)數(shù)，為飛艇動(dòng)導(dǎo)數(shù)的獲取提供一種新的思路。

1 動(dòng)導(dǎo)數(shù)CFD計(jì)算方法

1.1 縱向組合動(dòng)導(dǎo)數(shù)計(jì)算方法

飛艇俯仰組合動(dòng)導(dǎo)數(shù)是指俯仰力矩對(duì)俯仰角速度和俯仰力矩對(duì)迎角變化率的導(dǎo)數(shù)之和，可以通過模擬飛艇縱向俯仰振蕩獲取，假定飛行器繞重心俯仰振蕩方程為：

式中：θ為俯仰角，θ0為初始迎角，θm為振幅，ω為振蕩角頻率。

因此俯仰角速度可以寫為：

從飛艇運(yùn)動(dòng)狀態(tài)可知，飛艇迎角等于俯仰角，迎角變化率等于俯仰角速度，即：

飛艇俯仰力矩系數(shù)通過泰勒展開可以寫為：

帶入表達(dá)式（1）和（2）可以得到如下形式：

式中：Cm0為初始力矩，Cmα為縱向靜導(dǎo)數(shù)，Cmα˙為俯仰力矩對(duì)迎角變化率的導(dǎo)數(shù)，Cmq為俯仰力矩對(duì)俯仰角速度的導(dǎo)數(shù)，Cmα˙+Cmq為縱向組合動(dòng)導(dǎo)數(shù)。

從俯仰力矩的表達(dá)式可以看出，非定常模擬時(shí)間足夠長(zhǎng)時(shí)，俯仰力矩表現(xiàn)出和俯仰振蕩相同的周期性。

本文采用最小二乘法對(duì)縱向組合動(dòng)導(dǎo)數(shù)進(jìn)行參數(shù)辨識(shí)：根據(jù)縱向力矩的周期性，得到力矩周期性數(shù)據(jù)，并通過數(shù)據(jù)擬合得到如下表達(dá)式：

由此，通過擬合曲線可以辨識(shí)出縱向組合動(dòng)導(dǎo)數(shù)為：

升力系數(shù)組合動(dòng)導(dǎo)數(shù)可采用類似方法求解。

1.2 縱向單獨(dú)動(dòng)導(dǎo)數(shù)計(jì)算方法

飛艇縱向單獨(dú)動(dòng)導(dǎo)數(shù)包括俯仰力矩對(duì)俯仰角速度的導(dǎo)數(shù)Cmq和俯仰力矩對(duì)迎角變化率的導(dǎo)數(shù)Cmα˙。對(duì)于Cmq，可以通過旋轉(zhuǎn)流場(chǎng)法獲取，結(jié)合式（9）獲取的組合動(dòng)導(dǎo)數(shù)，進(jìn)而可以得到Cmα˙。

旋轉(zhuǎn)流場(chǎng)法即模擬飛艇的定常拉升運(yùn)動(dòng)，定常拉升運(yùn)動(dòng)是指飛艇在鉛錘面內(nèi)做角速度不變的圓周運(yùn)動(dòng)，由于速度始終沿著圓周的切線方向，故運(yùn)動(dòng)中迎角保持不變，俯仰角不斷變化。如果選擇固連在飛艇上的坐標(biāo)系，相對(duì)該坐標(biāo)系流動(dòng)就是定常的，因此可以進(jìn)行定常計(jì)算。

如圖1所示，保證飛艇速度不變，給定飛艇不同的運(yùn)動(dòng)角速度，可得到不同的俯仰力矩系數(shù)，俯仰阻尼導(dǎo)數(shù)表達(dá)式為

圖1 飛艇定常拉升運(yùn)動(dòng)

1.3 數(shù)值計(jì)算方法

計(jì)算采用非結(jié)構(gòu)瞬態(tài)滑移網(wǎng)格，滑移網(wǎng)格整個(gè)計(jì)算域分為靜止區(qū)域和滑移區(qū)域，兩個(gè)計(jì)算區(qū)域通過交界面interface進(jìn)行數(shù)據(jù)傳遞，滑移區(qū)域整體相對(duì)靜止區(qū)域進(jìn)行旋轉(zhuǎn)，運(yùn)動(dòng)過程中計(jì)算網(wǎng)格質(zhì)量不發(fā)生改變，非常適合飛行器的非定常氣動(dòng)計(jì)算。

控制方程為雷諾平均NS方程，湍流模型選用Realize k-e模型，k-e模型適合旋轉(zhuǎn)流動(dòng)、強(qiáng)逆壓梯度的邊界層流動(dòng)、流動(dòng)分離等情況，算法為simple算法，壓力項(xiàng)、對(duì)流項(xiàng)采用二階迎風(fēng)格式離散，擴(kuò)散項(xiàng)采用中心差分格式離散。

2 計(jì)算結(jié)果分析

2.1 Finner標(biāo)模驗(yàn)證

2.1.1 計(jì)算模型與計(jì)算條件

計(jì)算采用國際標(biāo)模Finner導(dǎo)彈模型，模型尺寸如圖2所示，模型重心到頭部距離為5d，d為彈體直徑。

圖2 Finner導(dǎo)彈模型

計(jì)算采用非結(jié)構(gòu)六面體笛卡爾網(wǎng)格，y+在1左右，網(wǎng)格質(zhì)量良好。采用非定常計(jì)算，計(jì)算Ma取風(fēng)洞試驗(yàn)馬赫數(shù)1.58，減縮頻率k取0.05，振蕩幅值取1°，初始迎角為0°。

2.1.2 計(jì)算結(jié)果分析

圖3是迎角和俯仰力矩系數(shù)一個(gè)周期內(nèi)隨時(shí)間變化曲線，可以看出迎角和俯仰力矩系數(shù)存在相位差，具有遲滯效應(yīng)。圖4為俯仰力矩系數(shù)遲滯環(huán)計(jì)算結(jié)果與文獻(xiàn)[5]計(jì)算結(jié)果的比較，可以看出，計(jì)算遲滯環(huán)曲線與文獻(xiàn)[5]基本一致。

圖3 俯仰力矩系數(shù)與迎角隨時(shí)間變化曲線

圖4 俯仰力矩系數(shù)遲滯環(huán)曲線

表1為本文計(jì)算結(jié)果與風(fēng)洞試驗(yàn)結(jié)果的對(duì)比，可以看出縱向和橫航向動(dòng)導(dǎo)數(shù)計(jì)算結(jié)果與試驗(yàn)結(jié)果誤差在5%以內(nèi)，表明本文計(jì)算方法可靠，動(dòng)導(dǎo)數(shù)計(jì)算精度較高。

表1 縱向動(dòng)導(dǎo)數(shù)計(jì)算結(jié)果

2.2 飛艇動(dòng)導(dǎo)數(shù)計(jì)算

2.2.1 計(jì)算模型與網(wǎng)格

計(jì)算對(duì)象為某型飛艇，計(jì)算模型和計(jì)算網(wǎng)格如圖5，計(jì)算網(wǎng)格采用非結(jié)構(gòu)六面體笛卡爾網(wǎng)格，網(wǎng)格質(zhì)量良好。靜止區(qū)域和滑移區(qū)域網(wǎng)格總數(shù)為430萬，計(jì)算模型周圍流場(chǎng)以及尾流區(qū)均進(jìn)行了不同程度的空間加密處理，不僅保證網(wǎng)格尺寸空間增長(zhǎng)的均勻過渡，又能對(duì)尾流進(jìn)行較好的捕捉。

圖5 計(jì)算模型和網(wǎng)格

2.2.2 計(jì)算條件

計(jì)算條件為標(biāo)準(zhǔn)海平面大氣條件，計(jì)算初始迎角為0°，速度20m/s，參考長(zhǎng)度50.176m，參考面積260.3m2，采用非定常計(jì)算。飛艇縱向低頻小振幅運(yùn)動(dòng)規(guī)律為：，振蕩頻率為1Hz，振蕩中心為飛艇重心，飛艇的俯仰振蕩通過編寫UDF實(shí)現(xiàn)。

2.2.3 計(jì)算結(jié)果及討論

圖6為采用俯仰振蕩法計(jì)算飛艇一個(gè)周期內(nèi)俯仰力矩系數(shù)隨計(jì)算迎角的變化曲線，即遲滯環(huán)。遲滯環(huán)為逆時(shí)針即非定常氣動(dòng)力做負(fù)功，縱向組合動(dòng)導(dǎo)數(shù)為負(fù)值。

圖6 遲滯環(huán)

給定飛艇無量綱化角速度分別為0.2189、0.4378、0.6568，圖7、圖8分別為俯仰力矩系數(shù)和升力系數(shù)隨俯仰角速度變化曲線。可以看出，俯仰力矩系數(shù)和升力系數(shù)隨俯仰角速度增加線性變化，通過線性擬合即可得到Cmq和CLq。

圖7 俯仰力矩系數(shù)隨俯仰角速度變化

圖8 升力系數(shù)隨俯仰角速度變化

表2為0度迎角下，飛艇縱向組合動(dòng)導(dǎo)數(shù)和單獨(dú)動(dòng)導(dǎo)數(shù)數(shù)值計(jì)算結(jié)果。

表2 縱向動(dòng)導(dǎo)數(shù)計(jì)算結(jié)果

由此可見，通過俯仰振蕩法和旋轉(zhuǎn)流場(chǎng)法可以獲得飛艇縱向組合動(dòng)導(dǎo)數(shù)和單獨(dú)動(dòng)導(dǎo)數(shù)，橫航向動(dòng)導(dǎo)數(shù)也可通過類似的方法獲取。由于飛艇和飛機(jī)氣動(dòng)外形差別較大，現(xiàn)有的飛機(jī)動(dòng)導(dǎo)數(shù)工程估算法并不適用飛艇，而俯仰振蕩法和旋轉(zhuǎn)流場(chǎng)法是飛艇動(dòng)導(dǎo)數(shù)求解的有效方法，具有較高的工程應(yīng)用價(jià)值。

3 結(jié)論

本文利用俯仰振蕩法對(duì)國際標(biāo)模Finner導(dǎo)彈進(jìn)行了動(dòng)導(dǎo)數(shù)計(jì)算，并將計(jì)算結(jié)果與風(fēng)洞試驗(yàn)結(jié)果進(jìn)行了對(duì)比。利用俯仰振蕩法和旋轉(zhuǎn)流場(chǎng)法對(duì)某型飛艇縱向組合動(dòng)導(dǎo)數(shù)和單獨(dú)動(dòng)導(dǎo)數(shù)進(jìn)行了計(jì)算，得到主要結(jié)論如下：①通過俯仰振蕩法對(duì)Finner導(dǎo)彈動(dòng)導(dǎo)數(shù)進(jìn)行了數(shù)值計(jì)算，計(jì)算結(jié)果與風(fēng)洞試驗(yàn)結(jié)果進(jìn)行了對(duì)比，誤差在5%以內(nèi)，說明本文數(shù)值方法的正確性。②通過俯仰振蕩法和旋轉(zhuǎn)流場(chǎng)法可以獲得飛艇縱向組合動(dòng)導(dǎo)數(shù)和單獨(dú)動(dòng)導(dǎo)數(shù)，為飛艇動(dòng)導(dǎo)數(shù)的獲取提供一種新的思路，具有較高的工程應(yīng)用價(jià)值。