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

（1. 航空工業(yè)成都飛機工業(yè)（集團）有限責任公司，成都610092；2. 中國氣象局氣象探測中心，北京100081）

1 引 言

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

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

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

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

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

2 物理模型

2.1 模擬對象

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

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

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

2.2 非定常RANS方程

采用雷諾平均N-S 方程求解流體域，非定常可壓縮RANS方程在直角坐標系下積分形式如下：

圖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)為切應力張量，包含分子粘性應力和雷諾應力。雷諾應力項需要通過湍流模式來封閉，在模擬中采用可實現(xiàn)的k-ε湍流模式[12]。

2.3 六自由度動力學方程

探空儀投放后在氣動力和重力的作用下自由運動，其運動遵循六自由度動力學方程，該方程在探空儀固連坐標系O1-X1Y1Z1下的表達式為[13]：

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

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

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

3 數(shù)值模擬方法

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

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

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

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

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

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

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

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

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

3.2 其他

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

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

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

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

4.1 非定常模擬收斂性驗證

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

子迭代步數(shù)分別取為10 步和15 步，模擬所得探空儀軌跡及姿態(tài)角分別如圖6、圖7 所示，其中位移以探空儀長度L無量綱化處理?？梢?，兩個模擬所得的探空儀軌跡吻合很好，歐拉角差別也很小，尤其是偏航角和俯仰角吻合很好，而滾轉(zhuǎn)角對探空儀氣動力影響小，因此為了兼顧計算效率，在模擬中選擇子迭代步數(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ǔ)上，選取時間步長分別為0.002s、0.001s和0.0005s，模擬所得探空儀軌跡及姿態(tài)角分別如圖8、圖9 所示?？梢?，不同時間步長模擬所得的探空儀軌跡吻合很好，而姿態(tài)角存在一定的差別。當時間步長大于0.001s時，姿態(tài)角的差別較為明顯，而時間步長小于0.001s 時，姿態(tài)角的差別較小，因此在一般模擬中選擇時間步長為0.001s。

圖8 不同時間步長模擬的探空儀軌跡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 探空儀運動分析

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

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

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

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

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

4.3 吊艙尾跡影響

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

圖12 探空儀所受氣動力Fig.12 The aerodynamic force of the dropsonde

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

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

4.4 流體粘性的影響

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

圖15 不同時刻吊艙尾跡流場Fig.15 The wake flow of the pod

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

圖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é) 論

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

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

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

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

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

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