陶善聰,呂 達(dá),張石玉,呂 蒙,時(shí)曉天,周江華

(1.中國航天空氣動(dòng)力技術(shù)研究院,北京 100074;2.中國科學(xué)院 光電研究院,北京 100094)

平流層飛艇在空中預(yù)警、通訊中繼、導(dǎo)航定位和環(huán)境監(jiān)測等方面有重要的軍用、民用價(jià)值,并具有應(yīng)用潛力,受到了研究人員的普遍關(guān)注。大型飛艇在平流層低速飛行時(shí),由于黏性效應(yīng)帶動(dòng)周圍的空氣一起運(yùn)動(dòng),由此產(chǎn)生空氣對飛艇的反作用力,可形成附加在飛艇表面的附加質(zhì)量或者附加力矩。附加質(zhì)量的獲取方法有工程估算、風(fēng)洞試驗(yàn)和數(shù)值模擬等手段,其中數(shù)值模擬因具有計(jì)算精度高、模擬成本小、周期短等優(yōu)勢而獲得廣泛應(yīng)用。

已有的研究結(jié)果表明,對于大型飛艇,其自身慣性所引起的附加質(zhì)量和力矩可能嚴(yán)重影響飛艇氣動(dòng)特性和操控特性。劉丹等通過Hess-Smith方法數(shù)值計(jì)算了飛艇在平流層運(yùn)動(dòng)時(shí)的附加質(zhì)量,指出附加質(zhì)量對飛艇運(yùn)動(dòng)具有較大影響,不能隨意假設(shè)。馬燁等對比分析了CFD(computational fluid dynamic,CFD)方法、理論方法和Hess-Smith方法計(jì)算標(biāo)準(zhǔn)模型橢球的附加質(zhì)量,結(jié)果表明CFD方法計(jì)算方便,精確度較高,并且將此方法推廣到計(jì)算復(fù)雜飛艇外形的附加質(zhì)量?？滤嚱艿冉柚こ坦浪愫蛯?shí)驗(yàn)辨識兩種方法分析平流層飛艇的俯仰附加轉(zhuǎn)動(dòng)慣量,結(jié)果表明兩者計(jì)算結(jié)果都可以滿足工程設(shè)計(jì)的要求,但經(jīng)驗(yàn)公式估算的結(jié)果相對偏小。張文琦等利用CFD方法考慮附加質(zhì)量作用下浮升混合飛艇的運(yùn)動(dòng)性能,指出附加質(zhì)量效應(yīng)對浮升混合飛行器的操縱性存在顯著影響。

本文基于六自由度非定常勻加速和重疊網(wǎng)格CFD技術(shù),首先針對橢球標(biāo)模開展了附加質(zhì)量計(jì)算方法驗(yàn)證。針對倒Y型尾翼布局的大型飛艇,開展了附加質(zhì)量和附加力矩的數(shù)值模擬,得到了飛艇平動(dòng)和轉(zhuǎn)動(dòng)的附加質(zhì)量和力矩;分析了布局形式對于沿軸、軸平動(dòng)的附加質(zhì)量、轉(zhuǎn)動(dòng)的附加力矩的影響差異;通過飛行動(dòng)力學(xué)仿真手段,考察了附加質(zhì)量和力矩對于飛艇開環(huán)控制響應(yīng)的影響。

1 數(shù)值方法及算例驗(yàn)證

1.1 數(shù)值方法

針對三維不可壓縮流動(dòng)的RANS(Reynolds-averaged Navier-Stokes,RANS)控制方程,本文采用有限體積法開展SST-湍流模型的非定常數(shù)值模擬,基于二階隱式方法實(shí)現(xiàn)時(shí)間推進(jìn)。采用背景網(wǎng)格靜止、內(nèi)域網(wǎng)格運(yùn)動(dòng)的重疊網(wǎng)格技術(shù)模擬飛艇的平動(dòng)和轉(zhuǎn)動(dòng),開展附加質(zhì)量和力矩的數(shù)值計(jì)算。附加質(zhì)量和力矩的計(jì)算公式為

(1)

(2)

式中:和分別為物體勻加速運(yùn)動(dòng)下的力與力矩,和分別為物體勻速運(yùn)動(dòng)下的力與力矩,和分別為平動(dòng)和轉(zhuǎn)動(dòng)的加速度。

1.2 算例驗(yàn)證

針對橢球標(biāo)模(2=20 m,2=2=5.333 m,=1.225 kg/m)開展附加質(zhì)量和力矩的數(shù)值計(jì)算,并通過與附加質(zhì)量的理論值和基于Hess-Smith方法獲得的預(yù)測結(jié)果進(jìn)行對比,驗(yàn)證計(jì)算方法的預(yù)測精度。計(jì)算網(wǎng)格如圖1所示,網(wǎng)格規(guī)模為200萬網(wǎng)格單元。

圖1 橢球網(wǎng)格示意圖

對于平動(dòng),開展來流速度1 m/s、0攻角的定常流場計(jì)算,在此基礎(chǔ)上賦予橢球一個(gè)1 m/s的勻加速度前向運(yùn)動(dòng),由式(1)獲得附加質(zhì)量。對于轉(zhuǎn)動(dòng),設(shè)置橢球在靜止流場中以1 rad/s勻加速轉(zhuǎn)動(dòng),轉(zhuǎn)動(dòng)中心為橢球中心,由收斂的力矩和式(2)求得其附加力矩。圖2為模型表面的壓力和對稱面上的流線分布圖。圖2(a)的加速前向運(yùn)動(dòng)結(jié)果表明流動(dòng)無分離,高壓區(qū)和低壓區(qū)分別出現(xiàn)在頭部和尾部。圖2(b)表明流場關(guān)于橢球中心成旋轉(zhuǎn)對稱分布,線速度與半徑成正比。高壓區(qū)出現(xiàn)在橢球的左下角和右上角,低壓區(qū)出現(xiàn)在左上角和右下角。

圖2 橢球勻加速運(yùn)動(dòng)

表1給出了本文計(jì)算值與理論值等結(jié)果的對比。由表1可見,當(dāng)前計(jì)算結(jié)果,平動(dòng)附加質(zhì)量與理論值誤差不超過0.3%,轉(zhuǎn)動(dòng)附加力矩預(yù)測誤差不超過4%,當(dāng)前的數(shù)值計(jì)算方法具有較高的可靠性。

表1 橢球附加質(zhì)量計(jì)算結(jié)果對比

2 飛艇附加質(zhì)量計(jì)算

2.1 飛艇模型及計(jì)算條件

圖3所示為本文所考察的倒Y型尾翼布局的大型飛艇示意圖。計(jì)算網(wǎng)格如圖4所示。圖4(a)為背景網(wǎng)格,域長、寬、高均為飛艇軸長的20倍,網(wǎng)格量為630萬,對應(yīng)靜止域的計(jì)算;圖4(b)為內(nèi)域網(wǎng)格,域長、寬、高均為飛艇軸長的4倍,網(wǎng)格量為520萬,對應(yīng)運(yùn)動(dòng)域的計(jì)算,靜動(dòng)域界面網(wǎng)格數(shù)據(jù)交換由插值方法得到。計(jì)算參數(shù)見表2,表中,為來流靜壓,為來流密度,為來流速度,為參考長度,為參考面積。

圖3 飛艇模型示意圖

圖4 飛艇重疊網(wǎng)格示意圖

表2 計(jì)算參數(shù)

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

圖5給出了不同攻角、側(cè)滑角條件下靜態(tài)計(jì)算結(jié)果,其中上面3幅圖分別為=0、=0、=0截面及艇身相對壓力等值線分布,下面3幅圖均為艇身=0和=0截線上相對壓力分布。對于攻角=0°和側(cè)滑角=0°的工況,表面高壓區(qū)主要集中在頭部和尾部,而低壓區(qū)集中在艇身過渡區(qū)和尾翼前緣,上下表面壓力分布基本一致。而對于有攻角和有側(cè)滑角工況,高壓區(qū)向迎風(fēng)面聚集,低壓區(qū)向背風(fēng)面聚集。3種流場中,尾翼表面均存在明顯的壓力變化。

圖5 靜態(tài)流場示意圖

附加質(zhì)量計(jì)算中,攻角設(shè)置=0°,首先考察了非定常時(shí)間步長對附加質(zhì)量計(jì)算結(jié)果的影響。設(shè)置時(shí)間步長Δ分別為0.01 s,0.001 s,0.000 1 s,對比分析附加質(zhì)量的預(yù)測結(jié)果。為了避免勻加速平動(dòng)或者轉(zhuǎn)動(dòng)時(shí)間過長引起非物理現(xiàn)象對計(jì)算結(jié)果的干擾,本文附加質(zhì)量計(jì)算值均為第50物理時(shí)間步的計(jì)算結(jié)果。平動(dòng)的勻加速計(jì)算在靜態(tài)均勻來流的基礎(chǔ)上開展,通過氣動(dòng)力差和式(1)得到附加質(zhì)量;勻加速轉(zhuǎn)動(dòng)計(jì)算中則直接開展非定常計(jì)算,靜態(tài)力矩作零處理,通過式(2)得到附加力矩。

圖6(a)表明,勻加速運(yùn)動(dòng)基本流動(dòng)狀態(tài)與靜態(tài)流場類似,流動(dòng)附著在飛艇表面并且過渡光滑平穩(wěn)。圖6(b)表明,飛艇在得到沿軸的加速度之后,表面高壓區(qū)向下表面轉(zhuǎn)移,低壓區(qū)向上表面轉(zhuǎn)移。圖6(c)表明,與軸勻加速運(yùn)動(dòng)類似,高壓區(qū)轉(zhuǎn)移至下表面,低壓區(qū)轉(zhuǎn)移至上表面,流體能夠附著在艇身表面流動(dòng),但尾翼表面流動(dòng)現(xiàn)象稍有差別。圖6(d)表明,由于三組尾翼在-平面對稱分布,因此繞軸勻加速旋轉(zhuǎn)時(shí),尾翼表面流場也較為規(guī)整光滑,高壓區(qū)附著在尾翼左側(cè)表面,低壓區(qū)附著在尾翼右側(cè)表面,流體粒子從高壓區(qū)出發(fā),勻加速旋轉(zhuǎn)運(yùn)動(dòng)至低壓區(qū)。圖6(e)表明,高壓區(qū)聚集在飛艇頭部的上表面和尾部的上表面,低壓區(qū)聚集在飛艇頭部的下表面和尾部的上表面,流體粒子從高壓區(qū)出發(fā)旋轉(zhuǎn)運(yùn)動(dòng)至低壓區(qū),此物理現(xiàn)象與轉(zhuǎn)動(dòng)方向相對應(yīng)。圖6(f)表明,由于艇身的對稱性,其表面流場與繞軸轉(zhuǎn)動(dòng)相似,而尾翼布置在和方向并不對稱,表面壓力分布存在差別。

圖6 飛艇勻加速運(yùn)動(dòng)流場示意圖

不同時(shí)間步長下飛艇勻加速運(yùn)動(dòng)氣動(dòng)參數(shù)收斂曲線如圖7所示。

圖7 不同時(shí)間步長下飛艇勻加速運(yùn)動(dòng)氣動(dòng)參數(shù)收斂曲線

沿、、軸勻加速平動(dòng)時(shí)不同時(shí)間步長下力系數(shù)收斂曲線分別如圖7(a)、圖7(b)、圖7(c)所示。圖中,為阻力系數(shù),為升力系數(shù),為側(cè)力系數(shù)。由圖可知力系數(shù)在幾個(gè)時(shí)間步內(nèi)即可收斂。對比穩(wěn)定后的計(jì)算結(jié)果,確定采用時(shí)間步Δ=0.001 s的結(jié)果用于附加質(zhì)量的計(jì)算,分別由式(1)得到沿軸平動(dòng)的=350.49 kg,沿軸平動(dòng)的=3 433.97 kg,沿軸平動(dòng)的=3 434.05 kg。

圖7(d)～圖7(f)分別給出了不同時(shí)間步長下繞軸、軸、軸勻加速旋轉(zhuǎn)時(shí)力矩系數(shù)收斂曲線。圖中,為滾轉(zhuǎn)力矩系數(shù),為偏航力矩系數(shù),為俯仰力矩系數(shù)。對比計(jì)算結(jié)果,選擇Δ=0.001 s的計(jì)算結(jié)果用于附加力矩的數(shù)值計(jì)算,分別由式(2)得到繞軸的附加力矩=23 133.49 kg·m,繞軸的附加力矩=1 701 351.90 kg·m,繞軸的附加力矩=1 700 477.71 kg·m。

對于附加質(zhì)量交叉項(xiàng),需要求解變參考點(diǎn)的力矩。壁面中心氣動(dòng)力對參考點(diǎn)的力矩為

(3)

式中:下標(biāo)“a”代表壓心坐標(biāo),下標(biāo)“mc”代表參考點(diǎn)坐標(biāo),參考點(diǎn)坐標(biāo)由勻加速運(yùn)動(dòng)方程給出。

圖8給出了沿軸勻加速平動(dòng)Δ=0.001 s的變體心俯仰力矩收斂曲線。從圖中看出,變參考點(diǎn)俯仰力矩基本能夠穩(wěn)定收斂,由式(2)可得=-11 929.66 kg·m。對于附加質(zhì)量分量,同樣采用上述方法求解。當(dāng)物理時(shí)間步長Δ=0.001 s時(shí)飛艇沿軸勻加速運(yùn)動(dòng)變參考點(diǎn)偏航力矩系數(shù)收斂曲線如圖9所示,從圖中可以看出,本次計(jì)算的偏航力矩系數(shù)在幾個(gè)物理時(shí)間步內(nèi)就快速收斂至固定值,由式(2)可得=12 053.71 kg·m。

圖8 俯仰力矩系數(shù)收斂曲線

圖9 偏航力矩系數(shù)收斂曲線

為比對附加質(zhì)量的計(jì)算結(jié)果,此處借助切片理論方法,利用旋轉(zhuǎn)橢球近似計(jì)算艇身附加質(zhì)量。采用近似橢球理論的估算方法得到飛艇的艇身部分附加質(zhì)量、、、、、,如表3所示,切片理論計(jì)算結(jié)果與本文數(shù)值計(jì)算方法所得到飛艇附加質(zhì)量中的艇身部分相對偏差小于5%,說明本文數(shù)值計(jì)算結(jié)果的合理性。

表3 附加質(zhì)量計(jì)算結(jié)果

3 飛艇開環(huán)動(dòng)力學(xué)響應(yīng)

本節(jié)利用前文的計(jì)算結(jié)果,分別建立有、無考慮附加質(zhì)量的飛艇動(dòng)力學(xué)模型,利用仿真分析手段研究附加質(zhì)量對飛艇開環(huán)動(dòng)力學(xué)響應(yīng)的影響。

針對本文的飛艇模型,做出如下假設(shè):該浮空器是剛體;該浮空器完全被理想氣體包圍;不考慮地球的曲率和自轉(zhuǎn)的影響。根據(jù)動(dòng)量和動(dòng)量矩定理,在浮空器的體坐標(biāo)系中建立的動(dòng)力學(xué)方程組:

(4)

式中:為動(dòng)量,為角速度,為動(dòng)量矩,為浮空器受到的全部外力在體坐標(biāo)系上的投影,為浮空器受到的全部外力矩在體坐標(biāo)系上的投影。其中,浮空器受到理想流體慣性力、黏性位置力、黏性阻尼力、浮力、重力、螺旋槳的推力和螺旋槳的失衡力矩的作用。將浮空器所受的全部外力和外力矩代入到式(4)中,則可得該浮空器考慮附加質(zhì)量的動(dòng)力學(xué)方程為

(5)

(6)

再由地面坐標(biāo)系與體坐標(biāo)系之間的轉(zhuǎn)換關(guān)系,得到浮空器的運(yùn)動(dòng)學(xué)方程為

(7)

將式中的附加質(zhì)量全部取為0,則可得到未考慮附加質(zhì)量的動(dòng)力學(xué)模型。建立該浮空器的考慮附加質(zhì)量的動(dòng)力學(xué)模型和未考慮附加質(zhì)量的動(dòng)力學(xué)模型,則可得到該浮空器針對不同控制輸入的響應(yīng)。仿真過程中設(shè)置的初始狀態(tài):高度為20 km,飛行速度為10 m/s,各類姿態(tài)角為0,且攻角側(cè)滑角為0,初始坐標(biāo)設(shè)為(0,0,0),仿真過程的基本設(shè)置:仿真時(shí)長為100 s,仿真步長為0.01 s,迭代方法為Runge-Kutta。圖10為高度通道的開環(huán)響應(yīng)對比,圖11為俯仰通道的開環(huán)響應(yīng)對比,圖12為偏航距離通道的開環(huán)響應(yīng)對比,圖13為偏航角通道的開環(huán)響應(yīng)對比。由圖可知,考慮附加質(zhì)量的浮空器動(dòng)力學(xué)模型的響應(yīng)速度相比未考慮附加質(zhì)量的浮空器的響應(yīng)速度更慢,由此可推斷,附加質(zhì)量能夠引起飛艇這類大體積浮空器的響應(yīng)遲滯,從而降低其操控性。

圖10 高度通道的開環(huán)響應(yīng)對比

圖11 俯仰通道的開環(huán)響應(yīng)對比

圖12 偏航通道的開環(huán)響應(yīng)對比

圖13 偏航角通道的開環(huán)響應(yīng)對比

4 結(jié)束語

本文采用六自由度勻加速運(yùn)動(dòng)和重疊動(dòng)網(wǎng)格CFD技術(shù)模擬了大型飛艇在平流層的勻加速平動(dòng)和轉(zhuǎn)動(dòng),得到了附加慣性質(zhì)量和附加慣性矩,分析了附加質(zhì)量對動(dòng)力學(xué)響應(yīng)的影響。得到以下結(jié)論:①附加質(zhì)量的CFD計(jì)算存在對時(shí)間步的依賴性,需要嚴(yán)格考察計(jì)算結(jié)果對于物理時(shí)間步收斂的無關(guān)性,提高附加質(zhì)量計(jì)算的可靠性。②尾翼為倒Y型特征的飛艇,沿軸和軸勻加速平動(dòng)的附加慣性質(zhì)量和繞軸、軸勻加速轉(zhuǎn)動(dòng)的附加慣性矩均近似相等。③在以往的飛行器設(shè)計(jì)當(dāng)中,附加質(zhì)量的影響通?？梢院雎?但是從動(dòng)力學(xué)響應(yīng)分析的角度說明,對于飛艇這類大體積、慢速度的飛行器來說,附加質(zhì)量引起的響應(yīng)遲滯現(xiàn)象不可忽略。