無人機(jī)下投探空儀投放過程數(shù)值模擬

肖良華，閆林明，曹世坤，郭 然，陳 斌

（1. 航空工業(yè)成都飛機(jī)工業(yè)（集團(tuán)）有限責(zé)任公司，成都610092；2. 中國(guó)氣象局氣象探測(cè)中心，北京100081）

1 引 言

利用飛機(jī)進(jìn)行下投探空，可對(duì)大氣溫度、壓力、濕度和風(fēng)速風(fēng)向等氣象要素的垂直分布進(jìn)行高分辨率測(cè)量，為臺(tái)風(fēng)的監(jiān)測(cè)和預(yù)報(bào)提供數(shù)據(jù)支持，因此具有極其重要的實(shí)際應(yīng)用價(jià)值。考慮到臺(tái)風(fēng)的特點(diǎn)及其探測(cè)的風(fēng)險(xiǎn)，高空長(zhǎng)航時(shí)無人機(jī)是下投探空的首選平臺(tái)[1-2]。為此，成都飛機(jī)工業(yè)（集團(tuán)）有限公司基于自主研制的云影無人機(jī)平臺(tái)設(shè)計(jì)了機(jī)載下投探空系統(tǒng)[3]。

在具體投放中，探空儀按照一定的時(shí)間間隔從無人機(jī)翼下吊艙尾部投出。類似“彈類”物體從飛機(jī)上的投放，探空儀的投放也需要開展安全性分析，以避免其投出后同飛機(jī)機(jī)體或尾翼發(fā)生碰撞。一般可以通過風(fēng)洞試驗(yàn)或者數(shù)值模擬的手段來研究探空儀投放過程。試驗(yàn)方法有動(dòng)力相似法、捕獲軌跡法和網(wǎng)格法[4]。由于探空儀重量輕、尺寸小，軌跡法和網(wǎng)格法試驗(yàn)難度大，而基于動(dòng)力相似法的縮比試驗(yàn)又難以保證弗勞德數(shù)（Fr）等參數(shù)的相似，如同時(shí)滿足重量、慣量及重心要求比較困難[5]，因此總體而言試驗(yàn)方法難度大、成本高。數(shù)值模擬則可以以較低成本、按照真實(shí)尺寸開展研究。

探空儀投放過程的數(shù)值模擬需要通過非定常流場(chǎng)的模擬來獲得氣動(dòng)力，然后將氣動(dòng)力代入動(dòng)力學(xué)方程中求解其速度及角速度，最后通過運(yùn)動(dòng)學(xué)方程積分獲得運(yùn)動(dòng)軌跡及姿態(tài)。顯然，氣動(dòng)力的求解尤為關(guān)鍵。投放過程中，探空儀同載機(jī)的相對(duì)位置一直在變，流體域的求解網(wǎng)格需要隨之變化。目前針對(duì)探空儀投放過程的數(shù)值模擬鮮見報(bào)道，其具體模擬方法可以參考“彈類”物體從飛機(jī)上分離的數(shù)值模擬。如張啟南等[6]基于面元法和工程估算法的氣動(dòng)力獲取方法開發(fā)了外掛物與載機(jī)分離過程的數(shù)值模擬系統(tǒng)。李孝偉等[7]利用動(dòng)態(tài)嵌套網(wǎng)格技術(shù)，通過求解Euler 方程和六自由度動(dòng)力學(xué)方程，模擬了副油箱從飛機(jī)上的投放過程。Deryl 等[8]、馮文梁等[9]采用基于Euler 方程的方法，模擬了武器外掛物的投放過程，模擬結(jié)果同風(fēng)洞試驗(yàn)結(jié)果吻合較好。楊磊[10]、杜小強(qiáng)[11]通過求解雷諾平均Navier-Stokes（RANS）方程獲得氣動(dòng)力，分別模擬了機(jī)載武器投放過程和其流動(dòng)特性。可見隨著計(jì)算能力的提升，氣動(dòng)力獲得方法越來越復(fù)雜，也越來越準(zhǔn)確。

同時(shí)，同“彈類”物體相比，探空儀的“密度”（重量同其外殼所占流體域體積之比）小得多，從而導(dǎo)致氣動(dòng)力在投放過程中的作用更加突出；轉(zhuǎn)動(dòng)慣量小、穩(wěn)定性較低，導(dǎo)致探空儀姿態(tài)更容易發(fā)生變化，因此對(duì)氣動(dòng)的模擬需要更加小心。而且探空儀為鈍頭體外形，開口端流動(dòng)分離嚴(yán)重，流場(chǎng)模擬難度大。

為此，本文采用動(dòng)態(tài)嵌套網(wǎng)格技術(shù)，通過非定常RANS 方程的求解獲得氣動(dòng)力，來實(shí)現(xiàn)探空儀投放過程的數(shù)值模擬。利用模擬結(jié)果分析了探空儀的運(yùn)動(dòng)及其受力變化過程。最后，為了評(píng)估流體粘性對(duì)模擬結(jié)果的影響，采用歐拉方法開展了同樣的數(shù)值模擬。

2 物理模型

2.1 模擬對(duì)象

模擬對(duì)象包括帶吊艙載機(jī)和探空儀。載機(jī)為云影無人機(jī)，其模型如圖1 所示，取坐標(biāo)系O-XYZ，坐標(biāo)原點(diǎn)為機(jī)頭位置，X軸為機(jī)體軸指向后，Z軸在對(duì)稱平面內(nèi)指向上，Y軸指向右，構(gòu)成右手坐標(biāo)系。模擬中，載機(jī)固定不動(dòng)，可選坐標(biāo)系O-XYZ為慣性坐標(biāo)系。

圖1 云影無人機(jī)模型及表面網(wǎng)格Fig.1 Model and surface mesh of Cloud Shadow unmanned aerial vehicle

探空儀模型如圖2所示，整體呈圓柱形，開口端內(nèi)安裝有傳感器（數(shù)值模擬的模型中忽略了內(nèi)部的傳感器），開口端的外支架可對(duì)傳感器起保護(hù)作用；封閉端內(nèi)安裝有降落傘等。投放時(shí)開口端先出艙，因此封閉端為迎風(fēng)面。為方便后續(xù)分析，建立如圖2 所示探空儀固連坐標(biāo)系O1-X1Y1Z1，原點(diǎn)O1為探空儀質(zhì)心，X1軸穿過探空儀軸線指向后，Y1軸及Z1軸與X1軸垂直，且構(gòu)成右手坐標(biāo)系。探空儀直徑D為65mm，總長(zhǎng)L為350mm，其中開口段長(zhǎng)24.8%L，質(zhì)量為380g，重心距離封閉端長(zhǎng)34.4%L。

2.2 非定常RANS方程

采用雷諾平均N-S 方程求解流體域，非定?？蓧嚎sRANS方程在直角坐標(biāo)系下積分形式如下：

圖2 探空儀模型Fig.2 Model of dropsonde

其中，ρ為空氣密度，p為壓力，u，v，w分別為氣流在x，y，z方向的速度分量，e為內(nèi)能，H1，H2，H3為無粘通量，J1，J2，J3為粘性通量，τij(i，j=x，y，z)為切應(yīng)力張量，包含分子粘性應(yīng)力和雷諾應(yīng)力。雷諾應(yīng)力項(xiàng)需要通過湍流模式來封閉，在模擬中采用可實(shí)現(xiàn)的k-ε湍流模式[12]。

2.3 六自由度動(dòng)力學(xué)方程

探空儀投放后在氣動(dòng)力和重力的作用下自由運(yùn)動(dòng)，其運(yùn)動(dòng)遵循六自由度動(dòng)力學(xué)方程，該方程在探空儀固連坐標(biāo)系O1-X1Y1Z1下的表達(dá)式為[13]：

其中，F(xiàn)→為探空儀受到的合力，包括重力和氣動(dòng)力；M→為探空儀所受相對(duì)質(zhì)心的氣動(dòng)力矩；V→、ω→分別為探空儀相對(duì)慣性坐標(biāo)系O-XYZ的速度和角速度，即相對(duì)于載機(jī)的速度和角速度；I為探空儀在質(zhì)心位置取的轉(zhuǎn)動(dòng)慣量，由于對(duì)稱性，慣性積Ixy=Ixz=Iyz=0，轉(zhuǎn)動(dòng)慣量

Ixx=1.93×10-4kg ?m3，Iyy=Izz=27.3×10-4kg ?m3。

通過式（1）迭代求解流體域，可以獲得探空儀所受氣動(dòng)力，將氣動(dòng)力帶入式（2），進(jìn)而可求得探空儀相對(duì)慣性坐標(biāo)系O-XYZ的速度和角速度，然后通過時(shí)間推進(jìn)求解運(yùn)動(dòng)學(xué)方程獲得下一時(shí)刻探空儀質(zhì)心位置及姿態(tài)角，其中姿態(tài)采用歐拉角偏航角Ψ、俯仰角Θ和滾轉(zhuǎn)角Φ表示，以更新后的相對(duì)位置重復(fù)上述步驟，即可獲得探空儀投放過程運(yùn)動(dòng)軌跡。

3 數(shù)值模擬方法

3.1 動(dòng)態(tài)嵌套網(wǎng)格方法

動(dòng)態(tài)嵌套網(wǎng)格方法具有網(wǎng)格生成簡(jiǎn)單、易于描述不同物體之間的相對(duì)運(yùn)動(dòng)的優(yōu)點(diǎn)，非常適合投放過程數(shù)值模擬[10]。

模擬中，分別生成載機(jī)網(wǎng)格（背景網(wǎng)格）和探空儀網(wǎng)格（子網(wǎng)格）。載機(jī)網(wǎng)格為靜止網(wǎng)格，其表面網(wǎng)格如圖1 所示，計(jì)算域取為20 倍機(jī)身長(zhǎng)度，非結(jié)構(gòu)網(wǎng)格總量約為一千萬，其中對(duì)探空儀運(yùn)動(dòng)過程中可能出現(xiàn)的區(qū)域進(jìn)行了加密處理。

探空儀網(wǎng)格為動(dòng)態(tài)網(wǎng)格，即探空儀網(wǎng)格相對(duì)探空儀靜止并隨探空儀一起運(yùn)動(dòng)，網(wǎng)格如圖3所示，網(wǎng)格量為110萬，在開口端尾跡區(qū)進(jìn)行了加密。

圖3 探空儀網(wǎng)格Fig.3 Mesh of dropsonde

一般情況下，使用探空儀網(wǎng)格外邊界切割載機(jī)網(wǎng)格，在切割邊界上進(jìn)行傳值，即可使整個(gè)流場(chǎng)域?qū)崿F(xiàn)非定常流場(chǎng)模擬。當(dāng)探空儀距離吊艙很近時(shí)（如圖4），探空儀網(wǎng)格同吊艙物面相交，此時(shí)需在切割載機(jī)網(wǎng)格前，先利用吊艙物面切割探空儀網(wǎng)格形成BC_cut0，如圖5紅色粗線所示，繼而利用BC_cut0和探空儀網(wǎng)格剩余外邊界一起切割載機(jī)網(wǎng)格，形成BC_cut1（圖4中藍(lán)色中等粗線）。

圖4 網(wǎng)格切割邊界Fig.4 Boundary of grid cutting

圖5 網(wǎng)格切割邊界局部放大圖Fig.5 Partial detail of grid cutting boundary

由于在模擬中的不同時(shí)刻，兩套網(wǎng)格的相對(duì)位置不同，因此網(wǎng)格的切割在每個(gè)時(shí)間步中都需要進(jìn)行。

3.2 其他

在CFD++軟件的數(shù)值模擬中除了采用動(dòng)態(tài)嵌套網(wǎng)格方法外，空間離散采用二階Total Variation Diminishing（TVD）格式，時(shí)間推進(jìn)則采用雙時(shí)間步隱式推進(jìn)方法，并采取多重網(wǎng)格加速措施。

邊界條件上，遠(yuǎn)場(chǎng)設(shè)置為基于特征值的進(jìn)出口邊界，來流馬赫數(shù)（Ma）為0.5，飛行高度設(shè)置為10km，攻角為-2°，側(cè)滑角為0°。物面邊界均設(shè)置為絕熱無滑移壁面。

實(shí)際投放中，探空儀以一定出射速度從吊艙中彈出。為方便仿真計(jì)算，將探空儀剛離開吊艙的狀態(tài)設(shè)為初始狀態(tài)，且假設(shè)此時(shí)探空儀軸線與吊艙軸線平行，出射速度沿探空儀軸線方向，大小為1m/s，初始時(shí)刻探空儀位置如圖4所示。由于探空儀出口軸線同X-Y平面成-1.2°夾角（迎風(fēng)端朝下），因此探空儀初始姿態(tài)歐拉角Θ=-1.2°。數(shù)值模擬時(shí)，先通過定常模擬獲得初始位置的穩(wěn)定流場(chǎng)，以便為自由運(yùn)動(dòng)模擬提供初始流場(chǎng)，然后開啟探空儀六自由度運(yùn)動(dòng)和非定常模擬。

4 仿真結(jié)果及分析

4.1 非定常模擬收斂性驗(yàn)證

非定常時(shí)間推進(jìn)過程中，隨著探空儀的運(yùn)動(dòng)，流場(chǎng)不斷變化，因此一方面需要確保流體域求解的收斂，另一方面需要時(shí)間推進(jìn)步長(zhǎng)足夠小，以保證探空儀運(yùn)動(dòng)軌跡積分過程的精確度。因此，為衡量非定常計(jì)算中的收斂情況，首先開展了不同子迭代步數(shù)和不同時(shí)間步長(zhǎng)的數(shù)值模擬。

子迭代步數(shù)分別取為10 步和15 步，模擬所得探空儀軌跡及姿態(tài)角分別如圖6、圖7 所示，其中位移以探空儀長(zhǎng)度L無量綱化處理。可見，兩個(gè)模擬所得的探空儀軌跡吻合很好，歐拉角差別也很小，尤其是偏航角和俯仰角吻合很好，而滾轉(zhuǎn)角對(duì)探空儀氣動(dòng)力影響小，因此為了兼顧計(jì)算效率，在模擬中選擇子迭代步數(shù)為10。

圖6 不同子迭代步數(shù)模擬的探空儀軌跡Fig.6 Dropsonde trajectory of simulations with different sub-iteration steps

圖7 不同子迭代步數(shù)模擬的探空儀歐拉角Fig.7 Dropsonde euler angle of simulations with different sub-iteration steps

在子迭代步數(shù)取為10的基礎(chǔ)上，選取時(shí)間步長(zhǎng)分別為0.002s、0.001s和0.0005s，模擬所得探空儀軌跡及姿態(tài)角分別如圖8、圖9 所示。可見，不同時(shí)間步長(zhǎng)模擬所得的探空儀軌跡吻合很好，而姿態(tài)角存在一定的差別。當(dāng)時(shí)間步長(zhǎng)大于0.001s時(shí)，姿態(tài)角的差別較為明顯，而時(shí)間步長(zhǎng)小于0.001s 時(shí)，姿態(tài)角的差別較小，因此在一般模擬中選擇時(shí)間步長(zhǎng)為0.001s。

圖8 不同時(shí)間步長(zhǎng)模擬的探空儀軌跡Fig.8 Dropsonde trajectory of simulations with different time step lengths

圖9 不同子迭代步數(shù)模擬的探空儀軌跡Fig.9 Dropsonde euler angle of simulations with different time step lengths

4.2 探空儀運(yùn)動(dòng)分析

基于子迭代步數(shù)為10、時(shí)間步長(zhǎng)為0.001s的數(shù)值模擬結(jié)果，作出不同時(shí)刻探空儀運(yùn)動(dòng)位置圖，如圖10 和圖11 所示。投放過程中，探空儀主要在流動(dòng)方向（向后）和重力方向（向下）運(yùn)動(dòng)，且流動(dòng)方向運(yùn)動(dòng)速度明顯大于重力方向。而一般的“彈類”物體的投放過程中，其運(yùn)動(dòng)主要體現(xiàn)在重力方向?？梢?，探空儀這類輕質(zhì)物體的投放過程受氣動(dòng)力的影響更為明顯。這也是本文選擇RANS方程而非歐拉方程求解流體域的原因（4.4 節(jié)中進(jìn)一步對(duì)比了RANS方法和歐拉方法模擬的結(jié)果）。

投放過程中，探空儀未出現(xiàn)翻滾，姿態(tài)角變化范圍不大，可見由于開口端的存在，使得探空儀氣動(dòng)力作用點(diǎn)處于重心之后，從而具有一定的穩(wěn)定性。探空儀離開吊艙0.2s 后則已經(jīng)脫離無人機(jī)機(jī)尾，因此數(shù)值模擬中只模擬0.25s時(shí)長(zhǎng)。

圖10 不同時(shí)刻探空儀位置（側(cè)視圖）Fig.10 Dropsonde release events(side view)

圖11 不同時(shí)刻探空儀位置（俯視圖）Fig.11 Dropsonde release events(top view)

進(jìn)一步對(duì)氣動(dòng)力進(jìn)行量化分析，作不同時(shí)刻三個(gè)方向的氣動(dòng)力和力矩如圖12、圖13所示?？梢?，y方向的氣動(dòng)力很小，其原因是偏航角Ψ較?。ㄈ鐖D7所示），并進(jìn)而使得y方向的速度（圖14）及位移均很小。俯仰角Θ較大，因此z方向氣動(dòng)力亦較大。至于x方向的氣動(dòng)力，即氣動(dòng)阻力，明顯大于其他兩個(gè)分量，其值可達(dá)重力（G）的5 倍以上，即x方向加速度可超過5g。對(duì)于氣動(dòng)力矩而言，由于探空儀呈圓柱體狀，因此其軸線方向（即x方向）氣動(dòng)力矩明顯小于其他兩個(gè)方向力矩。

4.3 吊艙尾跡影響

由圖12 可見，氣動(dòng)阻力在t=0.06s 之前的值較之后的值偏小，結(jié)合圖10 所示探空儀位置，可知探空儀在0.06s 之前距離吊艙較近，因此吊艙尾跡可影響探空儀受力。為此作出不同時(shí)刻過吊艙軸線的展向截面流場(chǎng)（吊艙尾跡附近，如圖15），可見t=0.06s 之前探空儀處于吊艙尾跡區(qū)內(nèi)，而尾跡區(qū)流速較低，故探空儀受到的氣動(dòng)力偏小；當(dāng)時(shí)間達(dá)到0.1s 后，探空儀已經(jīng)脫離吊艙尾跡，其受到的來流幾乎就是載機(jī)飛行來流，故在探空儀姿態(tài)變化不大、相對(duì)載機(jī)運(yùn)動(dòng)速度也不太大的時(shí)候，氣動(dòng)阻力穩(wěn)定在較大值附近。

圖12 探空儀所受氣動(dòng)力Fig.12 The aerodynamic force of the dropsonde

圖13 探空儀所受相對(duì)質(zhì)心的氣動(dòng)力矩Fig.13 The aerodynamic moment of the dropsonde

圖14 不同時(shí)刻探空儀速度Fig.14 The velocity of the dropsonde

4.4 流體粘性的影響

在“彈類”物體的投放過程模擬中，流體域采用歐拉（Euler）方程求解，獲得的投放物軌跡及姿態(tài)同實(shí)驗(yàn)值吻合較好[8-9]。為衡量流體粘性對(duì)探空儀軌跡預(yù)測(cè)的影響，本文還采用Euler 方程來代替RANS 方程求解流體域，模擬探空儀投放過程。Euler 方程本質(zhì)上是忽略了RANS 方程（1）中的粘性項(xiàng)。同時(shí)，壁面邊界改為滑移壁面，其他設(shè)置保持不變。

圖15 不同時(shí)刻吊艙尾跡流場(chǎng)Fig.15 The wake flow of the pod

探空儀軌跡如圖16所示，可見歐拉方法獲得的軌跡同RANS 方法差別很大，其本質(zhì)上是由兩種方法模擬的流場(chǎng)存在很大差異引起的，如圖17 所示，歐拉方法模擬的初始時(shí)刻流場(chǎng)同RANS方法模擬的流場(chǎng)存在很大差別。尤其是吊艙尾跡區(qū)，RANS 方法模擬的尾跡區(qū)較大，而歐拉方法由于沒有考慮粘性，其尾跡區(qū)范圍小得多，因此探空儀所受氣動(dòng)力亦差別很大。

圖16 歐拉方法及RANS模擬的探空儀軌跡Fig.16 Dropsonde trajectory of simulations with Euler method and RANS

圖17 歐拉方法模擬的吊艙尾跡Fig.17 The pod wake flow of the Euler simulation

5 結(jié) 論

本文針對(duì)下投探空儀從云影無人機(jī)翼下吊艙尾部投放安全性問題，采用非結(jié)構(gòu)動(dòng)態(tài)嵌套網(wǎng)格技術(shù)，通過求解RANS方程及六自由度動(dòng)力-運(yùn)動(dòng)學(xué)方程，數(shù)值模擬了探空儀的投放過程，分析了其運(yùn)動(dòng)軌跡、姿態(tài)變化及受力情況。作為對(duì)比，同時(shí)開展了基于歐拉方法的數(shù)值模擬。通過研究探索出了適合探空儀投放仿真的數(shù)值模擬方法，并得到了下面主要結(jié)論：

（1）對(duì)于嵌套網(wǎng)格切割問題，一般情況下使用子網(wǎng)格外邊界切割背景網(wǎng)格即可；但當(dāng)子網(wǎng)格同背景網(wǎng)格物面相交時(shí)，需要先采用背景網(wǎng)格物面對(duì)子網(wǎng)格進(jìn)行切割。

（2）探空儀的“密度”（重量同其外殼所占流體域體積之比）小，因此投放過程中氣動(dòng)力明顯大于重力，流向加速度可達(dá)5g以上。

（3）由于探空儀設(shè)計(jì)有開口端，使其氣動(dòng)力作用點(diǎn)處于重心之后，從而具有一定的穩(wěn)定性，投放過程中未出現(xiàn)翻滾，姿態(tài)角變化范圍不大，因此阻力方向氣動(dòng)力明顯大于另外兩個(gè)方向。

（4）相比于吊艙，探空儀尺寸較小，其投放初始段受吊艙尾跡影響很大，氣動(dòng)阻力明顯較小。

（5）同“彈類”物體相比，探空儀外形不具有流線型，且開口端流場(chǎng)復(fù)雜，如流體域采用Euler 方程模擬，則探空儀投放過程的儀軌跡偏離很大，表明忽略流體粘性不能獲得合理的氣動(dòng)力及力矩。在探空儀投放過程的模擬中必須采用RANS，甚至更精確的流動(dòng)模擬手段。