孫 博，高永衛(wèi)，魏斌斌

(1 西北工業(yè)大學(xué)航空學(xué)院，陜西 西安 710072；2 西安交通大學(xué)航天航空學(xué)院，陜西 西安 710049)

0 引言

無人機(jī)可以實(shí)現(xiàn)多機(jī)協(xié)同作戰(zhàn),執(zhí)行偵察、干擾和網(wǎng)絡(luò)中繼等任務(wù),經(jīng)濟(jì)可行性高,具有廣闊的軍事應(yīng)用前景[1-2]。由于無人機(jī)飛行速度和航時(shí)的限制,往往需要大型載機(jī)在戰(zhàn)區(qū)外對(duì)無人機(jī)實(shí)施投送,為戰(zhàn)區(qū)內(nèi)無人機(jī)協(xié)同作戰(zhàn)奠定基礎(chǔ)。載機(jī)對(duì)無人機(jī)的投送有以下特點(diǎn)：1)載機(jī)的動(dòng)能和勢能可以施加于無人機(jī),增加了無人機(jī)的初始速度和高度;2)投送過程中,無人機(jī)受載機(jī)影響強(qiáng)[3-4];3)相較于其他分離物,無人機(jī)的氣動(dòng)外形更復(fù)雜、升力和力矩特性的影響更顯著;4)無人機(jī)投送姿態(tài)、彈射速度等對(duì)整個(gè)分離過程有很大影響[5]。上述特點(diǎn)增加了無人機(jī)投送的技術(shù)難度。為保障無人機(jī)、載機(jī)的安全,必須研究無人機(jī)安全分離的條件。

在研究初期,往往希望能夠快速給出不同狀態(tài)下投放物離機(jī)軌跡,研究安全分離邊界,以確定方案是否成立[6]。該階段往往使用計(jì)算準(zhǔn)度滿足一定要求的簡便方法。這些簡便方法通常是準(zhǔn)定常方法[7],是把連續(xù)的非定常分離過程簡化成多個(gè)時(shí)間間隔很小的定常狀態(tài)進(jìn)行研究的方法。此類方法計(jì)算量少,盡管無法準(zhǔn)確計(jì)算出非定常效應(yīng),但是對(duì)于非定常效應(yīng)不顯著的情況可以快速給出投放物的分離軌跡。此類準(zhǔn)定?？焖俜椒ㄓ址Q為工程估算方法。

目前,典型的工程估算方法研究分為三類[8-19]：基于理論分析的、基于數(shù)值模擬(CFD)的和其他非典型方法?；诶碚摲治龅牡湫凸ぷ饔校?)范召林使用源匯模型和渦格模型模擬載機(jī)體積效應(yīng)和法向力效應(yīng),考慮載機(jī)在投放物上產(chǎn)生的干擾速度分布[8]。采用當(dāng)量迎角概念計(jì)算彈翼的氣動(dòng)力,但是該方法的一些假設(shè)和當(dāng)量迎角中干擾因子的計(jì)算尚有待進(jìn)一步驗(yàn)證[8]。與飛行試驗(yàn)相比,滾轉(zhuǎn)角的模擬大約有10°的偏差。2)美國國家航空航天局采用細(xì)長體理論計(jì)算外掛物的氣動(dòng)特性,建立了分離仿真系統(tǒng),給出了外掛物的分離軌跡[12]。由于細(xì)長體理論的限制,不能用于非細(xì)長體外掛物(如副油箱,無人機(jī)等)的分離軌跡研究。基于數(shù)值模擬的典型工作有：1)張啟南使用數(shù)值模擬方法計(jì)算載機(jī)的非均勻流場,采用面元法與工程估算相結(jié)合的計(jì)算方法計(jì)算載機(jī)非均勻流場中投放物的氣動(dòng)力。計(jì)算一條分離軌跡只需30 min,但是同樣存在對(duì)滾轉(zhuǎn)角模擬較差的情況[9],與風(fēng)洞CTS實(shí)驗(yàn)結(jié)果相比,投放物的滾轉(zhuǎn)角大約存在40°的偏差。2)AEDC(arnold engineering development center)使用CFD結(jié)合工程方法預(yù)測了跨音速下投放物與戰(zhàn)斗機(jī)的分離軌跡[13,17]。與飛行試驗(yàn)相比,在俯仰和偏航方向有較好的一致性,但是對(duì)滾轉(zhuǎn)方向的姿態(tài)角存在大約12°的偏差。其他非典型方法有：1)高曉光基于流場角度軌跡技術(shù)的思想,將風(fēng)洞試驗(yàn)得到的載機(jī)干擾流場制成數(shù)據(jù)庫,根據(jù)外掛物的坐標(biāo)值查詢載機(jī)干擾流場,得到相應(yīng)的流場系數(shù),進(jìn)而求得氣動(dòng)力和力矩系數(shù),最終獲得外掛物的運(yùn)動(dòng)狀態(tài);但只是將外掛物當(dāng)做質(zhì)點(diǎn)分析,定性的確定了干擾流場對(duì)外掛物姿態(tài)變化的影響,未做定量研究[10]。2)張冬云在均勻來流條件下計(jì)算了離機(jī)人員的阻力情況,并對(duì)人員的受力進(jìn)行了假設(shè)和簡化,建立了動(dòng)力學(xué)模型,提出了較為簡單的分離軌跡估算公式。但是將離機(jī)人員當(dāng)做質(zhì)點(diǎn),且不考慮載機(jī)干擾流場的影響,只能進(jìn)行簡單的定性分析[11]。綜上分析,目前已有的工程估算方法能夠在研究初期快速給出分離軌跡,但獲得的數(shù)據(jù)信息單一且準(zhǔn)度較低。文中認(rèn)為,分離軌跡估算準(zhǔn)度偏低的重要原因之一是滾轉(zhuǎn)角計(jì)算偏差較大。滾轉(zhuǎn)角估計(jì)不準(zhǔn)確將嚴(yán)重影響到投放物氣動(dòng)力估算的準(zhǔn)確性,從而導(dǎo)致分離軌跡的工程估算結(jié)果可信度較低。

另外,在研究方案階段,通常已分別獲取了大量載機(jī)流場和無人機(jī)的氣動(dòng)性能數(shù)據(jù),如果利用這部分?jǐn)?shù)據(jù)進(jìn)行軌跡預(yù)測,一方面不用增加額外的計(jì)算工作量,另一方面,又能夠改進(jìn)軌跡預(yù)測的準(zhǔn)確性,這將是改進(jìn)工程估算方法準(zhǔn)確性的一種新思路?；诖?文中發(fā)展了一種基于雙參考點(diǎn)的無人機(jī)分離軌跡高準(zhǔn)度快速預(yù)測方法。其中雙參考點(diǎn)取在無人機(jī)半模型的氣動(dòng)中心,考慮無人機(jī)左右流場的不均勻性帶來的滾轉(zhuǎn)力矩影響。經(jīng)驗(yàn)證,該方法利用現(xiàn)有數(shù)據(jù)能夠快速的計(jì)算無人機(jī)離機(jī)軌跡和姿態(tài),結(jié)果滿足工程應(yīng)用要求。

1 計(jì)算模型和計(jì)算方法

1.1 軌跡預(yù)測方法

1.1.1 載機(jī)和無人機(jī)參數(shù)

針對(duì)載機(jī)投放無人機(jī)的安全分離問題開展分離軌跡的快速預(yù)測方法研究。投放位置的大氣條件是海拔高度500 m、壓力95 461 Pa、密度1.167 3 kg/m3、重力加速度9.805 1 m/s2、溫度284.9 K;參考坐標(biāo)系為載機(jī)坐標(biāo)系O1X1Y1Z1,與載機(jī)固連,原點(diǎn)O1位于分離起始點(diǎn);O1X1平行于載機(jī)軸線,指向載機(jī)尾部;O1Y1在載機(jī)對(duì)稱平面內(nèi),垂直O(jiān)1X1向上;O1Z1垂直于載機(jī)對(duì)稱平面,指向左。無人機(jī)翼展0.3 m,外形為常規(guī)布局,其他參數(shù)見表1。

表1 無人機(jī)模型參數(shù)Table 1 UAV model parameters

載機(jī)翼展4.084 m,載機(jī)外形簡化為機(jī)翼加機(jī)身,研究中載機(jī)流場計(jì)算條件選取飛行馬赫數(shù)為0.45,載機(jī)的迎角為2°。無人機(jī)初始速度、初始位置和姿態(tài)(以下簡稱位姿)全部為0。

外掛物運(yùn)動(dòng)方程是包含3個(gè)旋轉(zhuǎn)運(yùn)動(dòng)和3個(gè)平移運(yùn)動(dòng)的六自由度方程。為了便于研究外掛物的分離過程,文中定義無人機(jī)坐標(biāo)系O2x2y2z2,與無人機(jī)固連,原點(diǎn)O2位于無人機(jī)質(zhì)心;O2x2平行于無人機(jī)軸線,指向無人機(jī)尾部;O2y2在無人機(jī)對(duì)稱平面內(nèi),垂直O(jiān)2x2向上;O2z2垂直于無人機(jī)對(duì)稱平面,指向左。文中在無人機(jī)坐標(biāo)系下求解運(yùn)動(dòng)方程。將無人機(jī)軸系的軌跡和姿態(tài)變化通過坐標(biāo)變換到載機(jī)軸系下表示。即分離軌跡結(jié)果均定義在載機(jī)軸系,可以直觀分析無人機(jī)與載機(jī)的分離安全性。

1.1.2 分離軌跡快速預(yù)測方法流程

文中方法包含數(shù)據(jù)準(zhǔn)備、無人機(jī)氣流角和動(dòng)壓的計(jì)算以及六自由度方程的求解。

1)載機(jī)流場準(zhǔn)備

基于求解RANS方程的CFD方法計(jì)算載機(jī)在給定來流條件下的干擾流場,計(jì)算結(jié)果如圖1所示。

圖1 載機(jī)干擾流場速度矢量圖(α=2°,β=0°,Ma=0.45,H=500 m)Fig.1 Velocity cloud diagram of carrier interference flow field (α=2°,β=0°,Ma=0.45,H=500 m)

2)無人機(jī)氣動(dòng)特性模型的建立

在已有的CFD結(jié)果基礎(chǔ)上,使用BP神經(jīng)網(wǎng)絡(luò)模型構(gòu)建無人機(jī)氣動(dòng)性能數(shù)據(jù)庫。樣本數(shù)據(jù)為無人機(jī)在干凈流場(無載機(jī))下的氣動(dòng)性能CFD結(jié)果,計(jì)算條件如表2所示。

表2 無人機(jī)CFD計(jì)算狀態(tài)參數(shù)變化Table 2 Changes in state parameters of UAV CFD calculation

由于載機(jī)流場中投放無人機(jī)的氣動(dòng)參數(shù)與迎角側(cè)滑角之間沒有精確的數(shù)學(xué)表達(dá)式,文中采用BP神經(jīng)網(wǎng)絡(luò)進(jìn)行無人機(jī)氣動(dòng)參數(shù)的插值。

根據(jù)神經(jīng)網(wǎng)絡(luò)相關(guān)文獻(xiàn)[19],隱藏層激活函數(shù)采Sigmoid函數(shù),輸出層激活函數(shù)采用Purelin函數(shù)。模型訓(xùn)練時(shí),隨機(jī)選擇70%用于模型訓(xùn)練,15%用于模型驗(yàn)證,15%用于模型測試。預(yù)測誤差如圖2所示。

圖2 BP神經(jīng)網(wǎng)絡(luò)模型預(yù)測誤差Fig.2 Prediction error of BP neural network model

預(yù)測誤差如圖2所示。在大多數(shù)的迎角、側(cè)滑角狀態(tài)下,輸出的氣動(dòng)力/力矩系數(shù)誤差約為5%。當(dāng)側(cè)滑角為0°時(shí),滾轉(zhuǎn)力矩系數(shù)和偏航力矩系數(shù)非常小,接近于0;當(dāng)迎角為0°時(shí),俯仰力矩系數(shù)也非常小。因此這兩類情況下的力矩系數(shù)相對(duì)誤差顯得稍大?？梢?文中訓(xùn)練的神經(jīng)網(wǎng)絡(luò)模型預(yù)測誤差小,滿足快速預(yù)測方法需求。

3)六自由度運(yùn)動(dòng)方程求解

使用四階Runge-Kutta方法[18]求解六自由度方程。

1.2 無人機(jī)參考點(diǎn)處氣流角與動(dòng)壓的計(jì)算方法

研究有兩個(gè)前提：1)無人機(jī)對(duì)載機(jī)流場的干擾量是小量。對(duì)于文中的投放情形而言,載機(jī)與無人機(jī)的展長之比為17∶1,可認(rèn)為無人機(jī)對(duì)載機(jī)流場的干擾小,無人機(jī)附近的流場主要取決于載機(jī)流場。2)載機(jī)干擾流場的非定常效應(yīng)不顯著。無人機(jī)初始分離位置處于載機(jī)腹部下方,相對(duì)于機(jī)翼/機(jī)身后方的流場,涉及的流場相對(duì)平穩(wěn)。無人機(jī)是以正常的姿態(tài)、零初始速度的自由投放。下面說明文中參考點(diǎn)位置的選取、氣流角和動(dòng)壓的計(jì)算方法。

1.2.1 參考點(diǎn)位置選取

文中的參考點(diǎn)沿展向?qū)ΨQ布置,如圖3所示(下文簡稱雙參考點(diǎn)),分別取考慮了無人機(jī)機(jī)翼和平尾等升力面影響的半模型氣動(dòng)中心。平均氣動(dòng)弦長的中點(diǎn)總是與機(jī)翼面積的中心相重合[20]。機(jī)翼的氣動(dòng)中心展向位置為0.075 m,平尾的氣動(dòng)中心展向位置為0.023 m。使用式(1)計(jì)算雙參考點(diǎn)的位置：

圖3 雙參考點(diǎn)位置示意圖Fig.3 Double reference point position diagram

L=k1L1+k2L2

(1)

式中：L1為機(jī)翼氣動(dòng)中心展向距離;L2為平尾氣動(dòng)中心展向距離;k1,k2為機(jī)翼和平尾對(duì)參考點(diǎn)展向位置的權(quán)重系數(shù)。理論上,升力面的升力決定了k1,k2的大小,為了便于計(jì)算,使用式(2)、式(3)對(duì)k1,k2進(jìn)行計(jì)算。

(2)

(3)

式(2)～式(3)中,S1和S2分別為機(jī)翼面積和平尾面積,S1=0.00764 m2,S2=0.00258 m2。經(jīng)計(jì)算：k1=0.75,k2=0.25,L=0.062 m。

雙參考點(diǎn)的x坐標(biāo)取無人機(jī)機(jī)翼平均氣動(dòng)弦長的四分之一處,z坐標(biāo)落在無人機(jī)坐標(biāo)系O2x2z2平面。應(yīng)該指出,文中的雙參考點(diǎn)不是傳統(tǒng)意義上的氣動(dòng)中心,而是經(jīng)過研究確定的適用于軌跡快速預(yù)測的特征位置。

1.2.2 氣流角和動(dòng)壓計(jì)算

分離軌跡快速預(yù)測方法的核心是載機(jī)干擾流場下無人機(jī)的氣動(dòng)力和力矩的確定,其關(guān)鍵是無人機(jī)氣流角和動(dòng)壓的計(jì)算。

為了研究方便,首先考慮采用代數(shù)平均求取載機(jī)干擾流場中的氣流速度V1的3個(gè)分量：

(4)

式中：n為參考域內(nèi)的網(wǎng)格節(jié)點(diǎn)數(shù);V1i為載機(jī)干擾流場中的氣流速度V1的i方向分量,i=1,2,3(1表示x方向,2表示y方向,3表示z方向,下同);v1i,j為第j個(gè)網(wǎng)格點(diǎn)上速度v1的i方向分量。

定義無人機(jī)運(yùn)動(dòng)速度為V2,V2i(i=1,2,3)分別為V2的3個(gè)分量。無人機(jī)質(zhì)心在載機(jī)干擾流場下的相對(duì)氣流速度V3的3個(gè)分量的計(jì)算公式為：

V3i=V1i-V2i

(5)

由V3的分速度及合速度求得質(zhì)心流動(dòng)角,根據(jù)質(zhì)心流動(dòng)角和無人機(jī)姿態(tài)角求得氣流角即無人機(jī)的迎角α、側(cè)滑角β,兩者可表示為：

α=θxy+θz
β=θxz+θy

(6)

動(dòng)壓q可表示為：

(7)

式中ρ為自由來流空氣密度。

偏差計(jì)算公式為：

(8)

參考域大小對(duì)氣流角、動(dòng)壓的影響如圖4所示。黑色的線表示參考域大小對(duì)代數(shù)平均求取等效值的影響。參考域球體半徑取值較小時(shí),參考域內(nèi)網(wǎng)格節(jié)點(diǎn)少,等效值波動(dòng)較大;參考域球體半徑較大時(shí)(R>0.15 m,無人機(jī)初始位置距離載機(jī)為0.15 m),參考域會(huì)與載機(jī)相交,導(dǎo)致等效值變化較大。當(dāng)參考域球體半徑在0.05～0.10 m時(shí),迎角、側(cè)滑角變化小于0.02°,速度大小變化小于0.01 m/s。因此,如果采用代數(shù)平均方式,參考域半徑取為0.08 m是合適的。

圖4 參考域大小對(duì)氣流角、動(dòng)壓的影響Fig.4 Influence of reference domain size on airflow angle and dynamic pressure

反比例權(quán)重和正比例權(quán)重方法對(duì)結(jié)果的影響如圖4和圖5所示。k=0為代數(shù)平均,k<0為反比例權(quán)重法,k>0為正比例權(quán)重法。由圖4可知,紅色的線更平穩(wěn),黃色的線更陡峭,但是不改變代數(shù)平均求取的穩(wěn)定區(qū)域大小。

圖5 不同加權(quán)方式對(duì)求取速度和氣流角的影響Fig.5 Influence of different weighting methods on velocity and airflow angle

在參考域半徑大小為0.08 m情況下,3種加權(quán)方式對(duì)求取速度大小和氣流角的影響如圖5所示。藍(lán)色、綠色、紅色的圓柱表示3種加權(quán)方式,對(duì)迎角的影響偏差最大不超過0.004°,對(duì)側(cè)滑角的影響偏差最大不超過0.02°,對(duì)速度的影響偏差最大不超過0.01 m/s。可見,不同加權(quán)方式對(duì)等效氣流角和速度的影響很小。因此,文中采用簡單代數(shù)平均計(jì)算參考域內(nèi)的等效氣流角和速度。

2 計(jì)算結(jié)果與討論

2.1 預(yù)測結(jié)果與CFD結(jié)果的對(duì)比

CFD計(jì)算采用URANS方法求解N-S方程得到給定時(shí)刻的分離物穩(wěn)態(tài)干擾氣動(dòng)力,結(jié)合剛體6DOF方程得到分離物下一時(shí)刻位置和姿態(tài)角,通過嵌套網(wǎng)格,使空間網(wǎng)格適應(yīng)分離物新的位置和姿態(tài)進(jìn)入下一步氣動(dòng)力計(jì)算循環(huán);最后,得到分離體的分離軌跡。通過該計(jì)算方法得到的計(jì)算結(jié)果,通過項(xiàng)目驗(yàn)收,與實(shí)際分離情況基本吻合,可用作工程估算方法的驗(yàn)證結(jié)果。CFD嵌套網(wǎng)格計(jì)算的無人機(jī)分離過程如圖6所示。

圖6 CFD嵌套網(wǎng)格計(jì)算的無人機(jī)分離過程Fig.6 UAV separation process based on CFD nested Grid computing

文中方法得出的無人機(jī)分離軌跡與CFD嵌套網(wǎng)格計(jì)算的分離軌跡對(duì)比如圖7所示。

圖7 快速預(yù)測方法與CFD計(jì)算結(jié)果對(duì)比Fig.7 Comparison between fast prediction method and CFD calculation result

由圖7可知,預(yù)測結(jié)果與CFD計(jì)算結(jié)果趨勢一致,準(zhǔn)確的得到了無人機(jī)分離過程的空間位置和姿態(tài)角變化。采用文中發(fā)展的預(yù)測方法計(jì)算的結(jié)果滿足工程要求。

2.2 參考點(diǎn)位置的影響

為了定量評(píng)價(jià)文中方法的預(yù)測偏差,進(jìn)一步驗(yàn)證文中求取參考點(diǎn)位置的合理性,文中將參考點(diǎn)位置范圍擴(kuò)大為由無人機(jī)質(zhì)心到翼尖,參考點(diǎn)在不同展向位置時(shí)軌跡預(yù)測的偏差如圖8所示。圖8(a)為空間位置偏差,圖8(b)為姿態(tài)角偏差。黑色實(shí)線、紅色實(shí)線和藍(lán)色實(shí)線分別表示x,y,z方向的位置和姿態(tài)?？招那虮硎举|(zhì)心參考點(diǎn)(單參考點(diǎn))、實(shí)心球表示使用2.2節(jié)方法計(jì)算的雙參考點(diǎn)。

圖8 不同展向位置的參考點(diǎn)預(yù)測分離軌跡姿態(tài)誤差Fig.8 Reference points at different spanwise positions predict separation trajectory posture error

圖8(a)中雙參考點(diǎn)較單參考點(diǎn)的x方向位置偏差增大了3.8×10-5m,y方向位置偏差減小了4.6×10-4m,z方向位置偏差增大了2.1×10-3m?？梢?參考點(diǎn)展向位置的變化對(duì)軌跡的空間位置偏差影響較小。圖8(b)中雙參考點(diǎn)較單參考點(diǎn)的滾轉(zhuǎn)角最大相對(duì)偏差從103%降低到3.4%,預(yù)測精度有了數(shù)量級(jí)上的提升,偏航角最大相對(duì)偏差從51.2%降低到16.5%,俯仰角最大相對(duì)偏差從43.1%降低到38.4%?？梢?雙參考點(diǎn)較單參考點(diǎn)能更好的預(yù)測無人機(jī)的姿態(tài)角,尤其是滾轉(zhuǎn)角。

文中方法較CFD方法耗時(shí)大大減少。在目前主流商用計(jì)算機(jī)上運(yùn)行該算法,從導(dǎo)入數(shù)據(jù)到獲得完整分離軌跡,文中的方法只需約120 s;而CFD方法給出一條完整分離軌跡需要3～5 d。當(dāng)計(jì)算幾十甚至上百個(gè)分離狀態(tài)時(shí),相較于CFD方法,該方法節(jié)省的時(shí)間將更可觀。

3 結(jié)論

發(fā)展了一種基于雙參考點(diǎn)的無人機(jī)離機(jī)軌跡和姿態(tài)高準(zhǔn)度快速預(yù)測方法,解決了傳統(tǒng)預(yù)測方法對(duì)滾轉(zhuǎn)角預(yù)測精度不高的不足,滿足研究初期快速獲得大量分離軌跡、確定安全邊界、制定后續(xù)實(shí)驗(yàn)與計(jì)算方案的需求。得出以下主要結(jié)論：

1)以無人機(jī)左右兩部分的氣動(dòng)中心為參考點(diǎn),考慮無人機(jī)左右機(jī)翼附近干擾流場的影響,計(jì)算由于流場的不對(duì)稱產(chǎn)生的滾轉(zhuǎn)氣動(dòng)力,對(duì)滾轉(zhuǎn)角的預(yù)測由傳統(tǒng)方法的103%降低到3.4%,預(yù)測精度有了數(shù)量級(jí)上的提升。

2)較以往取質(zhì)心處的流場信息發(fā)展的方法,文中以雙參考點(diǎn)為球心建立球形取值空間,通過區(qū)域內(nèi)多個(gè)數(shù)據(jù)的平均,增加了數(shù)據(jù)的平穩(wěn)性和準(zhǔn)確性,克服了采用單參考點(diǎn)數(shù)據(jù)隨機(jī)誤差大的缺陷。