陳宣合,張向強(qiáng),崔燕香

(1.中國科學(xué)院空天信息創(chuàng)新研究院,北京 100094;2.中國科學(xué)院大學(xué),北京 100049)

0 引言

對于浮空器的氣動(dòng)穩(wěn)定性,無動(dòng)力的系留氣球是靜穩(wěn)定的,中低空動(dòng)力飛艇是靜不穩(wěn)定的。臨近空間飛艇并非中低空飛艇的比例放大版,其顯著特點(diǎn)在于偏航控制機(jī)構(gòu)不同。中低空飛艇的偏航控制通過舵面提供控制力矩,該力矩與空速的平方成正比;而臨近空間的大氣密度僅為中低空的1/10,低速時(shí)舵面效率更低;此外,低空飛艇舵面機(jī)構(gòu)的重量、低溫低氣壓特性均不適合臨近空間飛艇。因此,臨近空間飛艇的偏航控制采用螺旋槳提供力矩。由于電機(jī)功率的限制,偏航或差動(dòng)螺旋槳的最大控制力矩基本固定,導(dǎo)致臨近空間飛艇空速超出設(shè)計(jì)值后,面臨控制力矩不足的局面,此時(shí),若臨近空間飛艇動(dòng)穩(wěn)定裕度較低,可能出現(xiàn)偏航通道發(fā)散,國內(nèi)近兩年的平流層飛行試驗(yàn)中已發(fā)生多次。高抗風(fēng)需求下的橫側(cè)向穩(wěn)定性問題需加強(qiáng)研究。

近年來,國外臨近空間飛艇飛行次數(shù)有所減少,臨近空間飛艇操控與氣動(dòng)穩(wěn)定性方面公開發(fā)表的研究成果很少。本文針對兩種已有的飛艇構(gòu)型,通過CFD仿真得出靜導(dǎo)數(shù)與動(dòng)導(dǎo)數(shù),進(jìn)行穩(wěn)定性分析,為分析偏航通道發(fā)散問題的機(jī)理提供依據(jù),同時(shí)為平流層飛艇氣動(dòng)設(shè)計(jì)提供參考。

1 模型與計(jì)算處理方法

1.1 飛艇模型

本文利用已有的HALE-D和HiSentinel兩種飛艇模型進(jìn)行數(shù)值仿真,模型示例如圖1所示,模型數(shù)據(jù)參看參考文獻(xiàn)[1]。使用ICEM CFD生成網(wǎng)格,附面層層數(shù)為30層,模型表面Y plus<1,空間六面體網(wǎng)格總數(shù)約800萬。

圖1 模型示例

坐標(biāo)系的定義不同于飛機(jī)等常規(guī)飛行器,艇身機(jī)體坐標(biāo)系原點(diǎn)為飛艇體心,方向定義:X向后沿機(jī)身縱軸;Z向左垂直于模型對稱面;Y向上垂直于XZ平面。計(jì)算輸出的氣流坐標(biāo)系原點(diǎn)為飛艇重心,方向定義:X向后沿遠(yuǎn)場氣流方向;Y向上垂直X軸并位于飛艇對稱面內(nèi);Z向左垂直于XY平面。

飛行速度相對于機(jī)體坐標(biāo)軸系的角度定義如下:側(cè)滑角為飛行速度與飛艇對稱面的夾角,當(dāng)飛行速度在橫軸分量為正時(shí)側(cè)滑角為正;迎角為飛行速度在飛艇對稱面的投影與縱軸的夾角,當(dāng)飛行速度沿豎軸分量為正時(shí)迎角為正。

側(cè)向力、所有力矩結(jié)果均為體軸系,升阻力、升阻比均為氣流坐標(biāo)系,動(dòng)導(dǎo)數(shù)結(jié)果均為體軸系。

機(jī)體坐標(biāo)系向氣流坐標(biāo)系轉(zhuǎn)換的公式為:

Oxayaza=By·Bz·Oxbybzb

(1)

其中:Oxayaza為氣流坐標(biāo)系坐標(biāo);Oxbybzb為機(jī)體坐標(biāo)系坐標(biāo);By為側(cè)滑角轉(zhuǎn)換矩陣;Bz為迎角轉(zhuǎn)換矩陣;α為迎角;β為側(cè)滑角。

1.2 穩(wěn)定性參數(shù)識(shí)別方法

本文通過Fluent軟件自編UDF程序,采用模型強(qiáng)迫振蕩方法來獲得動(dòng)導(dǎo)數(shù)[2-4]。定常流動(dòng)的NS方程采用SIMPLEC算法求解,非定常流動(dòng)的NS方程采用Coupled算法求解,計(jì)算模型采用SST-Kω模型。

(2)

其中:j為l或m,分別指升力和力矩系數(shù);V為來流速度。俯仰受迫振動(dòng)有以下運(yùn)動(dòng)學(xué)關(guān)系:

(3)

其中:αA為受迫周期性振蕩振幅;w為運(yùn)動(dòng)頻率。整理可得:

(4)

(5)

則

(6)

(7)

公式(6)為靜導(dǎo)數(shù),公式(7)為組合導(dǎo)數(shù),其它振動(dòng)方向處理方式與俯仰類似。

1.3 動(dòng)網(wǎng)格方法

模型強(qiáng)迫振蕩方法需要結(jié)合動(dòng)網(wǎng)格方法實(shí)現(xiàn),在CFD仿真軟件中,實(shí)現(xiàn)動(dòng)網(wǎng)格的常用方法有:剛性網(wǎng)格[5]、滑移網(wǎng)格[6]和重構(gòu)網(wǎng)格[7]。三種動(dòng)網(wǎng)格生成方式均會(huì)在計(jì)算域存在網(wǎng)格重疊,插值誤差會(huì)隨著計(jì)算時(shí)間推移而逐漸累積[8]。結(jié)合模型結(jié)構(gòu)網(wǎng)格劃分特點(diǎn),參考飛機(jī)動(dòng)導(dǎo)數(shù)及穩(wěn)定性相關(guān)文獻(xiàn),選用剛性網(wǎng)格結(jié)果,誤差相比另外兩種小[9]。

2 數(shù)值計(jì)算及分析

2.1 算例驗(yàn)證

為保證數(shù)值模擬結(jié)果可靠,首先,選用Finner[10-11]導(dǎo)彈模型驗(yàn)證穩(wěn)定性參數(shù)求解過程與精度。Finner導(dǎo)彈模型如圖2所示。

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

計(jì)算條件與參考文獻(xiàn)[12]一致,馬赫數(shù)Ma1.58,初始迎角α0=0°,減縮頻率k=0.015 822 6,振幅αm=1°,振蕩規(guī)律可寫為:

α=α0+αmsin(ωt)=1°sin(17t)

(8)

驗(yàn)證過程從初始迎角α0出發(fā)做定常計(jì)算,待殘差收斂及力矩系數(shù)穩(wěn)定在某值附近時(shí),開始做振蕩動(dòng)態(tài)計(jì)算。迎角及俯仰力矩系數(shù)隨時(shí)間變化圖如圖3所示,圖中顯示出了氣動(dòng)力矩的遲滯效應(yīng)。

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

根據(jù)積分法得到俯仰組合動(dòng)導(dǎo)數(shù),對照參考文獻(xiàn)[12]的風(fēng)洞試驗(yàn)結(jié)果,動(dòng)導(dǎo)數(shù)驗(yàn)證結(jié)果如表1所示。

表1 動(dòng)導(dǎo)數(shù)驗(yàn)證結(jié)果

結(jié)合參考文獻(xiàn)[12-13],本文驗(yàn)證相對誤差在3%以下,說明數(shù)值模擬結(jié)果比較準(zhǔn)確,方法可行。

2.2 縱向靜態(tài)氣動(dòng)力

圖4至圖8為迎角在-24°至24°范圍內(nèi),風(fēng)軸升力系數(shù)Cl、風(fēng)軸阻力系數(shù)Cd、升阻比K、體軸俯仰力矩系數(shù)Cmz和體軸法向力系數(shù)Cy的靜態(tài)參數(shù)變化圖。如圖4所示,HiSentinel(以下簡稱H1)的升力系數(shù)變化圖形似三次函數(shù),在0°迎角附近存在一段緩慢變化區(qū)域,且迎角繼續(xù)增大時(shí),升力系數(shù)也在均勻較快增長;HALE-D(以下簡稱H2)的升力系數(shù)變化圖在該范圍內(nèi)近似直線,升力系數(shù)增長變化明顯。

圖4 風(fēng)軸升力系數(shù)Cl靜態(tài)參數(shù)變化圖

在圖5所示的阻力系數(shù)變化圖中,H1與H2變化圖近似二次函數(shù),隨著迎角的增大,H2的阻力系數(shù)增長相較于H1更快。

圖5 風(fēng)軸阻力系數(shù)Cd靜態(tài)參數(shù)變化圖

圖6所示中,H1與H2升阻比變化圖形態(tài)基本相似,但H2的最大升阻比所對應(yīng)的迎角絕對值較H1略小。

圖6 升阻比K靜態(tài)參數(shù)變化圖

圖7所示中兩種模型俯仰力矩變化圖形似中心對稱正弦變化圖,在±24°內(nèi)都不穩(wěn)定,拐點(diǎn)約為-16°和16°。相較于H1模型而言,H2模型的力矩系數(shù)變化范圍較小。

圖7 體軸俯仰力矩系數(shù)Cmz靜態(tài)參數(shù)變化圖

由圖8可以看出,法向力系數(shù)變化圖趨勢與升力系數(shù)變化圖有較大相似性,相比于升力系數(shù)變化圖,法向力系數(shù)隨迎角變化更加明顯。

圖8 體軸法向力系數(shù)Cy靜態(tài)參數(shù)變化圖

2.3 橫向靜態(tài)氣動(dòng)力

圖9至圖11為側(cè)滑角在0°至24°范圍內(nèi),風(fēng)軸阻力系數(shù)Cd、體軸側(cè)向力系數(shù)Cz、體軸偏航力矩系數(shù)Cmy的靜態(tài)參數(shù)變化圖。

圖9 風(fēng)軸阻力系數(shù)Cd靜態(tài)參數(shù)變化圖

圖10 體軸側(cè)向力系數(shù)Cz靜態(tài)參數(shù)變化圖

圖11 體軸偏航力矩系數(shù)Cmy靜態(tài)參數(shù)變化圖

結(jié)合縱向與橫向的靜態(tài)穩(wěn)定性參數(shù)變化圖來看,H2模型相較于H1模型,隨相應(yīng)方向角度的增大,力的系數(shù)增長變化較為顯著,而在力矩系數(shù)上,兩種模型的變化趨勢大致相同,但H1模型的變化幅度更大。

2.4 俯仰動(dòng)穩(wěn)定性參數(shù)

在得到上述靜態(tài)穩(wěn)定性參數(shù)的基礎(chǔ)上,強(qiáng)迫飛艇振動(dòng),計(jì)算動(dòng)態(tài)穩(wěn)定性參數(shù)。兩種模型的強(qiáng)迫振動(dòng)頻率為0.1 Hz,振幅為3°,縱向俯仰力矩系數(shù)在0°迎角遲滯環(huán)曲線如圖12所示?？梢杂^測到,遲滯曲線以相應(yīng)的靜態(tài)穩(wěn)定性參數(shù)值為中心,呈橢圓形偏移。

圖12 縱向俯仰力矩系數(shù)在0°迎角遲滯環(huán)曲線

由遲滯環(huán)提取動(dòng)導(dǎo)數(shù),縱向振蕩俯仰力矩組合導(dǎo)數(shù)變化圖及縱向振蕩升力組合導(dǎo)數(shù)變化圖如圖13和圖14所示。整體上,模型俯仰振蕩俯仰力矩組合導(dǎo)數(shù)在整個(gè)迎角范圍內(nèi)為負(fù)值,模型為阻尼狀態(tài)。模型俯仰振蕩升力組合導(dǎo)數(shù)在整個(gè)迎角范圍內(nèi)為負(fù)值。

圖13 縱向振蕩俯仰力矩組合導(dǎo)數(shù)變化圖

圖14 縱向振蕩升力組合導(dǎo)數(shù)變化圖

2.5 偏航動(dòng)穩(wěn)定性參數(shù)[14]

偏航通道頻率與振幅和俯仰的相同,偏航力矩系數(shù)在0°迎角的遲滯環(huán)曲線如圖15所示。橫向振蕩偏航力矩組合導(dǎo)數(shù)變化圖和橫向振蕩偏航側(cè)向力組合導(dǎo)數(shù)變化圖分別如圖16、17所示。在縱向與橫向上,俯仰力矩組合導(dǎo)數(shù)及偏航力矩組合導(dǎo)數(shù)在對應(yīng)角度范圍內(nèi)均小于0,飛艇在對應(yīng)中心角度范圍附近為振幅收斂的振蕩,是動(dòng)穩(wěn)定狀態(tài)。

圖15 偏航力矩系數(shù)在0°迎角的遲滯環(huán)曲線

圖16 橫向振蕩偏航力矩組合導(dǎo)數(shù)變化圖

圖17 橫向振蕩偏航側(cè)向力組合導(dǎo)數(shù)變化圖

3 飛艇橫向穩(wěn)定性分析

3.1 飛艇橫向穩(wěn)定性與控制理論[1]

飛艇的橫向穩(wěn)定性是由飛艇的橫向運(yùn)動(dòng)和定向運(yùn)動(dòng)之間的相互作用決定的:橫向運(yùn)動(dòng)是由飛艇的側(cè)滑角產(chǎn)生的側(cè)向力造成的,當(dāng)飛艇橫向平移時(shí),側(cè)滑角在平移的相反方向上產(chǎn)生一個(gè)側(cè)力,側(cè)滑角的小擾動(dòng)會(huì)隨著時(shí)間的推移而減少;定向運(yùn)動(dòng)是由圍繞偏航軸旋轉(zhuǎn)產(chǎn)生的偏航力矩引起的。通常情況下,飛艇理論上是不穩(wěn)定的,因?yàn)榭諝鈩?dòng)力中心在參考點(diǎn)的前面,側(cè)滑角的小擾動(dòng)所產(chǎn)生的偏航力矩將傾向于增加該側(cè)滑角。側(cè)向和定向運(yùn)動(dòng)都受到飛艇側(cè)滑角和偏航率的影響,為了確定側(cè)向系統(tǒng)的整體穩(wěn)定性,必須同時(shí)考慮這兩個(gè)特性。

(9)

其中,V∞為遠(yuǎn)場速度。

側(cè)滑角速率簡化為體軸側(cè)向加速度除以真實(shí)空速。側(cè)向加速度等于側(cè)向力除以質(zhì)量加上由于體軸相對于慣性軸的旋轉(zhuǎn)而產(chǎn)生的科里奧利力。側(cè)向力用靜態(tài)側(cè)向力系數(shù)和動(dòng)態(tài)側(cè)向力偏航速率導(dǎo)數(shù)表示。

(10)

(11)

3.2 飛艇橫向穩(wěn)定性結(jié)果分析

穩(wěn)定性分析變化圖如圖18所示,Ⅰ區(qū)與Ⅱ區(qū)、Ⅱ區(qū)與Ⅲ區(qū)交界的兩條直線分別劃出H1和H2模型,理論上實(shí)現(xiàn)穩(wěn)定轉(zhuǎn)彎的穩(wěn)定區(qū)域和不穩(wěn)定區(qū)域。其中,Ⅰ區(qū)為H1的非穩(wěn)定轉(zhuǎn)彎區(qū)域,I區(qū)與Ⅱ區(qū)為H2模型轉(zhuǎn)彎的不穩(wěn)定區(qū)域。由于穩(wěn)定性和操縱性的孿生矛盾,H1模型即便存在不穩(wěn)定區(qū)域,相應(yīng)地也具備更好的操縱性;而H2模型即便距分界線較遠(yuǎn),穩(wěn)定裕度較大,相應(yīng)地飛行操縱性則大大降低。橫向偏航穩(wěn)定性需結(jié)合日后實(shí)際飛行試驗(yàn)進(jìn)一步加深研究。

圖18 穩(wěn)定性分析變化圖

4 結(jié)論

本文結(jié)合實(shí)際飛行需要,通過數(shù)值模擬方法得出兩種飛艇模型的靜態(tài)參數(shù)、動(dòng)穩(wěn)定參數(shù),并結(jié)合所得數(shù)據(jù)對飛艇偏航的穩(wěn)定性進(jìn)行分析,得出如下結(jié)論:

(1)靜態(tài)參數(shù)上,H1模型與H2模型的變化趨勢基本相同。在氣動(dòng)力系數(shù)上,H2模型隨迎角或側(cè)滑角的變化更加顯著;在氣動(dòng)力矩系數(shù)上,H1和H2在計(jì)算角度內(nèi)均靜不穩(wěn)定,H1模型變化幅值較H2模型的更大。此外,相同迎角(側(cè)滑角)時(shí),H2的阻力較大。

(2)動(dòng)穩(wěn)定性參數(shù)上,H1模型與H2模型變化范圍均較小,俯仰振蕩、偏航振蕩組合導(dǎo)數(shù)在整個(gè)角度變化范圍內(nèi)均為負(fù)值,模型處于阻尼狀態(tài)。

(3)橫向穩(wěn)定性分析表明,H2模型穩(wěn)定區(qū)間大于H1模型。

按照對常規(guī)飛行器運(yùn)動(dòng)穩(wěn)定性的理解,氣動(dòng)靜不穩(wěn)定即意味著運(yùn)動(dòng)不穩(wěn)定,只有依靠控制系統(tǒng)才可保持運(yùn)動(dòng)穩(wěn)定[15-16],但從已有的飛艇氣動(dòng)資料發(fā)現(xiàn),許多飛艇是靜不穩(wěn)定的[17]。相關(guān)文獻(xiàn)[18]指出,由于飛艇獨(dú)特的氣動(dòng)外形,動(dòng)導(dǎo)數(shù)對飛艇的氣動(dòng)影響是常規(guī)飛行器的數(shù)百倍,對飛艇穩(wěn)定性的影響極大。國內(nèi)平流層飛艇項(xiàng)目整體尚在小規(guī)模實(shí)驗(yàn)階段,技術(shù)方法有待完善。HALE-D和HiSentinel系列是國外飛艇經(jīng)過飛行實(shí)驗(yàn)的典型構(gòu)型,也是倒“Y”型尾翼和“X”型尾翼的典型代表。對于飛艇橫向穩(wěn)定性的研究,是以尾翼形狀為主體的,艇體外形的影響微乎其微。

由于穩(wěn)定性與操控性的孿生矛盾,飛艇的設(shè)計(jì)需根據(jù)實(shí)際需求,綜合考慮阻力、穩(wěn)定裕度和操縱性能之間的平衡。本文所研究的兩種不同尾翼形狀,可以根據(jù)不同的任務(wù)需求,應(yīng)用于飛艇控制系統(tǒng)設(shè)計(jì)中。