高速下集束薄片初始分離過程仿真研究

王政偉，王 浩，阮文俊，寧惠君，王金龍，王鵬新

（1.南京理工大學(xué)能源與動力工程學(xué)院，江蘇南京 210094；2.中國白城兵器試驗中心，吉林白城 137001）

高速下集束薄片初始分離過程仿真研究

王政偉1，＊，王 浩1，阮文俊1，寧惠君1，王金龍1，王鵬新2

（1.南京理工大學(xué)能源與動力工程學(xué)院，江蘇南京 210094；2.中國白城兵器試驗中心，吉林白城 137001）

為了研究諸如美國MJU-B50系列面源式紅外干擾彈圓柱形集束薄片在高速運動下發(fā)射的分離散布過程，建立了ALE有限體積法描述下的三維可壓縮RANS方程，耦合剛體運動學(xué)方程、碰撞檢測及響應(yīng)模型，數(shù)值模擬了載機馬赫數(shù)0.8條件下集束薄片特征段的三維非定常初始分離過程。研究結(jié)果表明：集束薄片分離初期碰撞頻繁，擾亂了薄片的一致性運動，加速了薄片分離過程，但是薄片的分離次序主要受薄片間的多體干擾影響，驗證了動態(tài)拋撒試驗過程中出現(xiàn)的薄片次序分離特性。

兵器科學(xué)與技術(shù)；多體分離；集束薄片；可壓縮流；流固耦合

0 引 言

隨著新型可成像式雷達(dá)的陸續(xù)裝備，傳統(tǒng)的點源式干擾彈越來越滿足不了作戰(zhàn)要求，面源式紅外干擾彈作為一種新型的紅外干擾系統(tǒng)應(yīng)運而生。如美國MJU-50B紅外干擾彈［1］，它采用自燃材料，從發(fā)射筒到點火前完全封閉，未氧化金屬薄片從圓筒中彈出后與空氣接觸逐漸氧化并輻射熱量。薄片的散布規(guī)律直接影響了干擾彈的作戰(zhàn)效果。

目前研究薄片云團的運動多采用建立誘餌彈整體質(zhì)點運動模型，認(rèn)為誘餌薄片為空間均勻球狀分布［2］或者采用正態(tài)隨機數(shù)方法［3］獲得云團分布。陳乃光［4］建立了低速下薄片的質(zhì)點模型，忽略了薄片分離過程中相互間氣動干擾，得出了自燃薄片的分離散布規(guī)律，但是薄片尺寸較大，受氣動作用明顯，并不適合顆粒質(zhì)點模型。黃蓓［4］建立了ALE有限體積法下的三維不可壓縮方程，耦合剛體運動方程數(shù)值模擬了薄片的低速下的靜態(tài)拋撒過程，并與試驗結(jié)果做了對比，其研究都僅限于靜態(tài)低速下的拋撒過程。但是在實戰(zhàn)中，干擾彈都是在戰(zhàn)機飛行過程中發(fā)射。為研究這種干擾彈在動態(tài)拋撒過程中薄片初始分離的形態(tài)，并解釋試驗中發(fā)現(xiàn)的薄片初始分離過程中頭尾次序分離的現(xiàn)象，本文基于剛體動力學(xué)理論和碰撞檢測、響應(yīng)模型，耦合流體計算方法對高速運動下集束薄片的分離過程進行了仿真分析。

1 數(shù)值模擬方法

1.1 薄片動態(tài)拋撒過程描述

圖1為面源式干擾彈工作示意圖。大數(shù)量圓柱狀薄片摞齊裝填，薄片在發(fā)射時成串射出，快速散開燃燒形成干擾云團。這個過程中薄片云團的散布形態(tài)與散布面積是干擾彈能否產(chǎn)生有效干擾的關(guān)鍵之一，所以研究薄片集束在戰(zhàn)機飛行過程中發(fā)射的薄片散開機理也是指導(dǎo)其改進設(shè)計的重要理論依據(jù)。

圖1 面源式干擾彈工作示意圖Fig.1 Structure diagram of surface-type infrared decoy

1.2 控制方程及離散方法

建立ALE有限體積法描述下的三維可壓縮方程［5］，其積分表達(dá)式為：

式中：Ω為控制體體積，S為控制體表面邊界，n為控制體邊界外法向單位矢量。守恒變量Q、分量E、F、G及Ev、Fv、Gv的表達(dá)式分別為：

式中：u、v、w為地面坐標(biāo)系下流體速度；ug、vg、wg為網(wǎng)格運動速度，wg；p、ρ、e為控制體內(nèi)的壓力、密度和單位質(zhì)量氣體的總能；qx、qy、qz為微元熱流量

湍流模型采用雙方程Realizable k-ε模型［5］，此模型對于不同雷諾數(shù)下的邊界層流動和帶有分離的流動等均有較好的適應(yīng)性。引入Reynolds平均法得到時均化的流動控制方程，將Reynolds應(yīng)力表示為湍流粘性系數(shù)的函數(shù)，湍流粘性系數(shù)由湍動能k和湍流耗散率ε求得。對于近壁區(qū)域低雷諾數(shù)流動情況，采用壁面函數(shù)法進行修正求解。

對控制方程及湍流方程用二階迎風(fēng)格式進行空間離散，時間離散采用隱式方案。薄片采用固壁邊界，外流場區(qū)域來流采用自由來流邊界，出流邊界由內(nèi)場外推處理，并根據(jù)湍動強度和特征長度來估算入口邊界處的k和ε值。

采用彈簧近似法和局部網(wǎng)格重構(gòu)相結(jié)合的非結(jié)構(gòu)動網(wǎng)格技術(shù)實現(xiàn)分離過程中物體邊界運動后流體計算網(wǎng)格的更新［6］。耦合剛體六自由度運動方程計算獲得邊界位移量［7－12］。

1.3 剛體運動耦合計算方法與碰撞判定方法

薄片在空氣中的運動過程中主要受氣動力作用和薄片間的接觸碰撞作用影響。本文碰撞判定定義為薄片之間的間距小于警戒距離時為碰撞作用。薄片間的碰撞定義為剛體的完全彈性碰撞，只有薄片的速度和角速度在碰撞后發(fā)生改變。

1.3.1 剛體碰撞檢測模型

薄片在分離初始位置時間距較小，分離過程中薄片的運動差異性會使其在分離過程中的碰撞頻率非常高。為有效減小數(shù)據(jù)計算量，提高檢測速度，把薄片間的碰撞判定分為兩個步驟：

（1）預(yù)判定。運用球體碰撞模型，以包裹薄片的最小球體為單元，若兩球體質(zhì)心間距離大于其半徑和，即，判定為薄片未碰撞；若兩球體質(zhì)心間距離小于等于其半徑和，即（r1＋r2），則進入第二步判定。

（2）精確判定。通過搜索兩薄片表面網(wǎng)格的質(zhì)心，尋找兩薄片間的最近距離，在聯(lián)系薄片間的相對速度、計算時間步長等參數(shù)建立起碰撞條件為：

式中，dij為t時刻兩薄片間的最小距離；u為t時刻兩薄片間的相對速度；Δt為此時的時間步長；k為修正系數(shù)，根據(jù)網(wǎng)格質(zhì)量等因素適當(dāng)設(shè)置；dmin為碰撞警戒距離。若t時刻滿足方程（2），則判定為碰撞，啟動碰撞響應(yīng)程序，若不滿足則判定為不碰撞跳到下一步流場計算。

1.3.2 薄片碰撞模型［6，13－14］

薄片的碰撞模型采用三維剛體碰撞模型，其基本假設(shè)如下：

（1）碰撞薄片均為剛體，碰撞前后薄片無變形，能量無損失；

（2）碰撞瞬間完成，薄片無相對滑動；

（3）碰撞瞬間，薄片間撞擊應(yīng)力遠(yuǎn)大于氣動力，不考慮碰撞時氣動力的影響。

（4）碰撞過程中不考慮多體間的碰撞，多體碰撞簡化為兩兩相碰。

碰撞模型如圖2所示。假設(shè)兩個薄片a和b，質(zhì)量分別為ma和mb，碰撞前a、b的質(zhì)心速度分別為va0和vb0，其體坐標(biāo)系下的角速度分別為ω、ω，其質(zhì)心到碰點的距離分別為ra、rb。坐標(biāo)系轉(zhuǎn)換張量為I，則a、b在碰撞前瞬間碰點的速度為：

圖2 碰撞模型示意圖Fig.2 Diagram of collision model

碰撞點的法向由薄片間的相對位置決定，碰撞后薄片分別受到一個沿法向的相等相反的沖量Im，作用在2個薄片上的沖量矩分別為Impn×ra和－Impn× rb，其中pn為碰撞點法向的單位矢量。根據(jù)動量和動量矩定理，物體碰撞時刻的平動和轉(zhuǎn)動方程分別為：

式中，va、vb為兩個剛體碰撞后的速度；ωa、ωb為薄片碰撞后的角速度；Ia、Ib為薄片轉(zhuǎn)動慣量。由彈性碰撞理論給出輔助方程：

式中，k為彈性恢復(fù)因數(shù)，對于完全彈性碰撞，k＝1；ua、ub分別為碰撞后2個薄片的速度。

由方程（3）～（9），可得出沖量表達(dá)式為：

聯(lián)合以上方程即可得薄片a、b碰撞后的速度和角速度。

1.3.3 薄片運動方程［15］

以固定坐標(biāo)系oxyz為地面坐標(biāo)系，以薄片質(zhì)心建立Cξηζ為薄片的彈體坐標(biāo)系，ψ、θ、φ為物體相對彈體平動坐標(biāo)系的歐拉角，則根據(jù)剛體運動方程可列出薄片質(zhì)心運動方程，用矩陣形式表述為：

式（10）中，m為薄片質(zhì)量矩陣，r為薄片質(zhì)心的坐標(biāo)矢量，r＝［x y z］T，F(xiàn)為薄片所受空氣動力在地面坐標(biāo)系下的坐標(biāo)矩陣為薄片所受空氣動力在薄片彈體坐標(biāo)系下的坐標(biāo)矩陣，G為重力，Λ為方向余弦陣。式（11）中，J為薄片的轉(zhuǎn)動慣量矩陣，ω為薄片角速度在彈體坐標(biāo)系下的坐標(biāo)矩陣，M 為薄片所受空氣動力矩在彈體坐標(biāo)系下的坐標(biāo)矩陣，

式（12）中q為歐拉角坐標(biāo)陣，q＝［ψ θ φ］T，K為與ω之間的轉(zhuǎn)換矩陣，

聯(lián)合方程（10）～（12）可知，如果已知薄片初速度、及所受空氣動力即可確定薄片的運動軌跡。

1.3.4 流固耦合方法及碰撞加入方法

采用弱耦合的方法分別獨立求解流體方程和運動方程，通過流固邊界交互數(shù)據(jù)。在一個計算步長內(nèi)首先求解流體方程。流體控制方程把計算出來的氣動力傳遞給固體運動方程，固體運動方程根據(jù)氣動力把計算出來的固體邊界位置和速度傳遞給流體控制方程和網(wǎng)格控制程序，網(wǎng)格控制程序控制固體邊界的移動、轉(zhuǎn)動并對畸變率大的網(wǎng)格進行重構(gòu)等等。碰撞檢測程序在網(wǎng)格控制程序之后運行。碰撞檢測程序經(jīng)過預(yù)判定后，若無法確定是否發(fā)生碰撞啟動精確判定程序。精確判定程序調(diào)用薄片邊界所有網(wǎng)格點位置變量，求解兩薄片間的最小間距，若最小間距小于警戒距離，則啟動碰撞模型，得出新的薄片運動參數(shù)并返回給流體方程與固體運動方程。否則，進行下一步計算。

2 數(shù)值計算與結(jié)果分析

2.1 計算模型與網(wǎng)格

在高速動態(tài)下的薄片分離過程中非定常作用明顯，薄片間的擾流復(fù)雜，氣動力作用明顯，碰撞頻繁，與低速靜態(tài)拋撒分離過程差異很大。

數(shù)值模擬的薄片模型取直徑d＝50mm，質(zhì)量為1g的均質(zhì)薄片，與拋撒試驗?zāi)Ｐ鸵恢拢穸葹?mm。薄片間距為2mm。如圖3所示，5片薄片在發(fā)射筒中以35m／s速度射出，受到飛機牽連速度為240m／s，根據(jù)相對速度原理，在流場仿真計算中給定流場速度x方向為－240m／s，y方向為－35m／s，薄片集束自下而上分別編號為1?！?＃。初始計算網(wǎng)格如圖4（a）所示，薄片間間隔采用加密處理，其最小尺寸為1／200d，其他部分網(wǎng)格用比例函數(shù)放大，為保證網(wǎng)格質(zhì)量，放大系數(shù)不宜超過1.2。

圖3 薄片拋撒模型示意圖Fig.3 Diagram of computational models

2.2 多薄片的分離過程與仿真分析

多薄片的分離流場仿真是包含多個運動邊界的雙向流固耦合問題，而且這種運動邊界是厚徑比很小的薄片狀物體，在高速流場中加速快，速度高，對網(wǎng)格劃分和網(wǎng)格重構(gòu)要求較高，必須以較小的時間步長來進行仿真。兼顧工作站計算能力和流場計算的準(zhǔn)確性，通過對小數(shù)量的薄片分離流場進行仿真，來對薄片的分離過程進行仿真研究。

圖4為五薄片分離過程xy截面的網(wǎng)格重構(gòu)變化圖。薄片附近網(wǎng)格采用函數(shù)加密處理，網(wǎng)格重構(gòu)均勻，計算模型初始網(wǎng)格為200萬，能很好地反映出流場的變化。隨著動邊界的移動和轉(zhuǎn)動，邊界層網(wǎng)格經(jīng)過彈簧近似與重構(gòu)，網(wǎng)格數(shù)量會逐漸增加，在2ms后達(dá)到800萬。更新后的網(wǎng)格能夠滿足流場計算條件，所獲得的流場計算結(jié)果收斂程度較好。

圖4 五片分離過程xy截面網(wǎng)格重構(gòu)變化圖Fig.4 Gird of five plate reconstruction course cross the xysection

圖5為五個薄片分離過程各時刻流場流線與壓力云圖，反映了多薄片在分離過程中位置變化和姿態(tài)。結(jié)合圖6薄片的速度變化曲線與圖7薄片運動的位移時間曲線可以看出，在xy方向上薄片多體的干擾明顯，脫離主體后減速過程迅速，其加速度大小和薄片自身所處位置和姿態(tài)關(guān)系最為密切。不受主體段遮擋氣流的脫離薄片，最大加速度可達(dá)100 000g以上。頭部薄片在氣動力的作用下從集束上方脫離，其余薄片從集束尾部脫離。從尾部分離的薄片分離次序規(guī)律，其位移變化趨勢相似，通過統(tǒng)計分析得出相鄰薄片間分離時間間隔約為0.6ms。在z軸方向上薄片的運動受碰撞作用的影響，速度改變量較大，但在集束薄片初始分離階段其位移與速度變化相對較小，所以在初始分離階段以xy方向的分離為主。

圖5 五片多體分離過程xy截面壓力云圖Fig.5 Pressure cloud imagery of five plate multibodies separation course cross the xysection

圖8為相鄰薄片最小間距時間。碰撞接觸雖然在分離初期1ms內(nèi)頻繁發(fā)生，結(jié)合圖6薄片運動速度曲線可以看出，碰撞產(chǎn)生的速度階躍在x、y運動方向并不明顯，并沒有大幅度改變整體運動趨勢，在薄片z運動方向上速度階躍較大，但是相對于x、y方向上的運動，其速度低、位移小，不足以改變其運動趨勢。所以碰撞并不是分離的主要原因，碰撞并沒有改變薄片分離的次序，但是碰撞產(chǎn)生的角速度與速度變化導(dǎo)致薄片間的差異性運動加快了薄片的分離過程。薄片多體間干擾明顯，脫離主體后減速過程迅速。尾部薄片分離次序規(guī)律基本遵循依次逐片分離的次序。

圖6 薄片運動速度時間圖Fig.6 Speed of plates change with time

圖7 薄片運動位移時間圖Fig.7 Displacement of plates change with time

圖8 相鄰薄片最小間距-時間圖Fig.8 Minimum distance between adjacent plate change with time

3 結(jié) 論

通過建立ALE有限體積法描述下的三維可壓縮方程，耦合剛體運動方程，并加入碰撞判定與碰撞響應(yīng)模型，數(shù)值模擬了高速下集束薄片的三維非定常初始分離過程。通過仿真分析得出以下結(jié)論：

（1）把碰撞問題加入到了流固耦合計算中，通過仿真模擬了高速下薄片群自干擾彈發(fā)射筒彈射后的三維非定常初始分離過程，該仿真模型可用于包含多個運動邊界的流固耦合與碰撞問題的流場問題求解。

（2）分析了薄片高速下多體分離中流場變化與薄片姿態(tài)的運動及薄片間的最小間距變化，薄片集束分離初期碰撞頻繁，擾亂了薄片的一致性運動，加速薄片分離過程，但是薄片的分離次序主要受薄片間的多體干擾影響。

（3）特征段多薄片分離次序基本穩(wěn)定，在空間上薄片集束頭部先脫離，尾部遵循依次脫落的分離過程，并通過統(tǒng)計分析得出了時間上的分離時序規(guī)律，可為薄片云團的仿真提供理論基礎(chǔ)。在下一步的工作中我們將根據(jù)仿真得出的薄片分離規(guī)律來對薄片云團的散布分離過程做進一步的仿真。

［1］ Gan Y L，Jiang C，Liu Y J，et al.Development of foreign airborne IR decoytechnology［J］.Electro-Optic Technology Application，2013，28（6）：13-17.（in Chinese）淦元柳，蔣沖，劉玉杰，等.國外機載紅外誘餌技術(shù)的發(fā)展［J］.光電技術(shù)應(yīng)用，2013，28（6）：13-17.

［2］ Fu X H，F(xiàn)an Q L.Research on modeling and simulation of special material decoy［J］.Electro-Optic Technology Application，2013，28（6）：81.（in Chinese）付曉紅，樊秋林.面源紅外誘餌的建模與仿真研究［J］.光電技術(shù)應(yīng)用，2013，28（6）：81.

［3］ Zhao F Y，Lu S，Jiang C，Research on the modeling method of surface-type infrared decoy［J］，Electro-Optic Technology Application，2012，27（2）：66.（in Chinese）趙非玉，盧山，蔣沖.面源紅外誘餌仿真建模方法研究［J］.光電技術(shù)用，2012，27（2）：66.

［4］Chen N G.Moving characters simulation of foils cloud used in optics-Electronic countermeasure［J］.Aerospace Electronic Warfare，2008，24（5）：20.（in Chinese）陳乃光.光電干擾箔片云團運動規(guī)律的仿真研究［J］.航天電子對抗，2008，24（5）：20.

［5］ Wang F J.The analysis of computational fluid dynamics［M］.Beijing：Tsinghua University Press，2011.（in Chinese）王福軍.計算流體動力學(xué)分析［M］.北京：清華大學(xué)出版社，2011.

［6］ Huang B，Wang H，Tao R Y.Numerical study of multi-platesseparating flow and characteristics［J］.Acta Aerodynamica Sinica，2013，31（2）：213-218.（in Chinese）黃蓓，王浩，陶如意.薄片分離過程流場特性的數(shù)值仿真研究［J］.空氣動力學(xué)學(xué)報，2013，31（2）：213-218.

［7］ Zhu B，Zhu X P，Zhou Z，et al.Simulation of unsteady multibody flowfield involving relative movement based on unstructured mesh［J］.Acta Aerodynamica Sinica，2013，31（2）：181-185.（in Chinese）朱冰，祝小平，周洲，等.基于非結(jié)構(gòu)動網(wǎng)格的多體分離數(shù)值仿真研究［J］.空氣動力學(xué)學(xué)報，2013，31（2）：181-185.

［8］ Duan X P，Chang X H，Zhang W P.A CFD and 6-DOF coupled solver for multiple moving object problems based on dynamic hybrid gird［J］.Acta Aerodynamica Sinica，2011，29（4）：447-452.（in Chinese）段旭鵬，常興華，張來平.基于動態(tài)混合網(wǎng)格的多體分離數(shù)值模擬方法［J］.空氣動力學(xué)學(xué)報，2011，29（4）：447-452.

［9］ Batina J T.Implicit flux-split Euler schemes for unsteady aerodynamic analysis involvingun structured dynamic meshes［R］.NASA TM-102732，1990.

［10］Guo Z.Numerical method with multi-body forunsteady flows with moving boundaries［D］.National University of Defense Technology，2002.（in Chinese）郭正.包含運動邊界的多體非定常流場數(shù)值模擬方法研究［D］.國防科技大學(xué)，2002.

［11］Zhang Y D，Ji C Q.Numerical methods with dynamic structured meshes forunsteady flows with moving boundaries［J］.Chinese Jouranal of Computational Physics，2006，23（2）：165-170.（in Chinese）張玉東，紀(jì)楚群.包含動邊界的非定常流場動網(wǎng)格數(shù)值模擬［J］.計算物理，2006，23（2）：165-170.

［12］Banerjee A K，Dickens J M.Dynamics of an arbitrary flexible body in large rotation and translation［J］.Journal of Guidance，Control and Dynamics，1990，13（2）：221-227.

［13］Huang B，Wang H，Wang S，et al.Flow field characteristics of multi-plates interference in descent［J］.Journal of Ballistics，2012，24（1）：41-46.（in Chinese）黃蓓，王浩，王帥，等.薄片狀體沉降過程中的多體干擾流場特性［J］.彈道學(xué)報，2012，24（1）：41-46.

［14］Lin J Z.Multiphase fluid dynamics of extraordinary particle［M］.Beijing：Science Press，2008.（in Chinese）林建忠.超常顆粒多相流體動力學(xué)［M］.北京：科學(xué)出版社，2008.

［15］Pu F.Exterior ballistics［M］.National Defense Industry Press，1980.（in Chinese）浦發(fā).外彈道學(xué)［M］.國防工業(yè)出版社，1980.

Simulation of plates group initial separation in high speed

Wang Zhengwei1，＊，Wang Hao1，Ruan Wenjun1，Ning Huijun1，Wang Jinglong1，Wang Pengxin2
（1.Nanjing University of Science and Technology，Nanjing 210094，China；2.China Baicheng Weapon Test Centre，Baicheng 137001，China）

In order to study the separation press of the cylindrical paltes launched form the surface type IR decoy such as MJU-B50armed in U.S.air force，the method solving threedimensional compressible RANS equations described by ALE finite volume method coupled with the rigid motion equations was established，and the collision detection and the collision response model were developed to solve the collision between the plates.The three-dimensional unsteady initial separation process of the several plates was numerical simulated at Mach number of 0.8.Research shows that the plates collide frequently at incipient motion under close initial locations，cause the inconsistent motion and accelerate the separation course.The separation order of the plates group is mainly affected by the aerodynamics of multi-plates interference，the plates separate in a fixed law and the verified order separation characteristics existed in the dynamic test of plates dispersal.

ordnance science and technology；multi-body separation；compressible flow；plates group；fluid-solid coupling

V411.8；V448.25＋3

：Adoi：10.7638／kqdlxxb-2014.0076

2014-07-23；

2014-12-02

王政偉＊（1987-），男，湖北荊門，博士研究生，研究方向：兵器發(fā)射理論與技術(shù).E-mail：wangzhengwei2421＠163.com

王政偉，王浩，阮文俊，等.高速下集束薄片初始分離過程仿真研究［J］.空氣動力學(xué)學(xué)報，2015，33（6）：828-834.

10.7638／kqdlxxb-2014.0076 Wang Z W，Wang H，Ruan W J，et al.Simulation of plates group initial separation in high speed［J］.Acta Aerodynamica Sinica，2015，33（6）：828-834.

0258-1825（2015）06-0828-07