徐學(xué)文，肖支才，曲 凱

（海軍航空大學(xué)，山東 煙臺(tái) 264001）

0 引言

1 仿真計(jì)算方程

飛機(jī)發(fā)射/投放彈是現(xiàn)代化戰(zhàn)爭(zhēng)中重要的作戰(zhàn)活動(dòng)和軍事打擊手段。隨著第5代戰(zhàn)斗機(jī)的研制成功及軍事應(yīng)用，為提高戰(zhàn)斗機(jī)的隱身功效，戰(zhàn)斗機(jī)均采用內(nèi)埋艙貯式發(fā)射導(dǎo)彈與投放炸彈。

在飛機(jī)打擊目標(biāo)過程中，飛機(jī)艙彈分離過程非常復(fù)雜：在飛機(jī)彈艙處產(chǎn)生流動(dòng)分離、激波干擾等復(fù)雜的流動(dòng)現(xiàn)象，在艙外表面還存在強(qiáng)烈的氣流剪切層[1]；飛機(jī)與導(dǎo)彈/炸彈分離之后，導(dǎo)彈/炸彈將進(jìn)行六自由度運(yùn)動(dòng)，飛行姿態(tài)極易受外界氣流、作用力和激波的影響產(chǎn)生較大變化，這不僅對(duì)載機(jī)的安全性造成嚴(yán)重的威脅，而且還對(duì)投放安全、姿態(tài)穩(wěn)定與打擊精度造成影響[2]。因此，研究飛機(jī)艙彈分離過程，提高艙彈分離品質(zhì)與導(dǎo)彈/炸彈的打擊精度，具有重要的軍事意義。

當(dāng)前，國(guó)內(nèi)外對(duì)艙彈分離問題的研究主要采用風(fēng)洞試驗(yàn)[3-4]、飛行試驗(yàn)[5-6]和數(shù)值仿真[7-8]3 種方法，其中，風(fēng)洞試驗(yàn)和飛行試驗(yàn)存在研究成本高、威脅性大、耗時(shí)長(zhǎng)等缺點(diǎn)。

近年來(lái)，隨著計(jì)算機(jī)性能提高和動(dòng)網(wǎng)格仿真技術(shù)進(jìn)步[9-11]，數(shù)值仿真已成為科學(xué)研究的主要手段。因此，本文采用計(jì)算流體動(dòng)力學(xué)（CFD）和剛體動(dòng)力學(xué)（RBD）模型耦合求解的數(shù)值仿真方法[12]，基于動(dòng)網(wǎng)格技術(shù)模擬飛機(jī)艙彈分離過程，研究炸彈投放分離后俯仰角度的變化。

1.1 RBD方程

飛機(jī)在投放彈時(shí)保持勻速直線飛行，然后打開艙門，把要投放的炸彈從彈艙彈射出去。拋射瞬間，炸彈受到重力FG、拋射力FS和拋射力矩MS的作用，脫離貯存架，炸彈做六自由度運(yùn)動(dòng)。艙彈分離后，炸彈還要受到來(lái)流氣體壓力FP、黏性摩擦力Fu和氣動(dòng)力矩Md的作用。在炸彈運(yùn)動(dòng)的數(shù)值模擬過程中，不考慮炸彈的材料特性（假設(shè)為剛體），將炸彈在空間的運(yùn)動(dòng)看作是導(dǎo)彈質(zhì)心的移動(dòng)和導(dǎo)彈繞質(zhì)心的轉(zhuǎn)動(dòng)的合成[13]，在慣性坐標(biāo)系(x,y,z)下可用線速度vc(vx,vy,vz)和角速度ωc(ωx,ωy,ωz)來(lái)描述，且它們都是時(shí)間t的函數(shù)。同時(shí)，在炸彈上建立以質(zhì)心為原點(diǎn)的隨動(dòng)體坐標(biāo)系(xb,yb,zb)，如圖1所示，ocxb為彈體軸線，指向頭部，oczb在彈體中心對(duì)稱面內(nèi)垂直ocxb軸線，指向下方。

圖1 炸彈受力及坐標(biāo)系Fig.1 Force and coordinate system of bomb

慣性坐標(biāo)系下炸彈質(zhì)心平移運(yùn)動(dòng)方程為：

式（1）中：v?c為慣性坐標(biāo)系下炸彈質(zhì)心的加速度；m為炸彈質(zhì)量；F為炸彈質(zhì)心處所受外力。

體坐標(biāo)系下炸彈角運(yùn)動(dòng)方程為：

式（2）中：ω?b為體坐標(biāo)系下炸彈的角轉(zhuǎn)動(dòng)加速度；L為轉(zhuǎn)動(dòng)慣量；Mb為力矩矢量。體坐標(biāo)系炸彈角運(yùn)動(dòng)方程可以通過坐標(biāo)轉(zhuǎn)換矩陣I轉(zhuǎn)換成慣性坐標(biāo)系下運(yùn)動(dòng)方程。

在仿真計(jì)算中，通過對(duì)式（1）（2）積分就可以確定炸彈質(zhì)心的位置和運(yùn)動(dòng)方向。設(shè)和分別表示當(dāng)前第n時(shí)間步質(zhì)心的位置和方向，則下一個(gè)時(shí)間步(n+1)時(shí)質(zhì)心的位置及方向?yàn)椋?/p>

剛體的位置矢量根據(jù)瞬時(shí)角速度ωc轉(zhuǎn)動(dòng)來(lái)確定，對(duì)于有限的轉(zhuǎn)動(dòng)角Δθ= |ωc|?Δt，炸彈位置矢量xr相對(duì)于質(zhì)心xc表示為：

式（5）中，eθ、er為單位矢量。

1.2 CFD方程

本文涉及對(duì)RBD的數(shù)值模擬，因而使用任意拉格朗日-歐拉方法（Arbitrary Lagrangian-Eulerian，ALE）描述的N-S 方程，對(duì)流場(chǎng)進(jìn)行描述[14-16]。在CFD 慣性坐標(biāo)系下，對(duì)于邊界移動(dòng)的任意控制體積V上的標(biāo)量φ（質(zhì)量ρ、速度u、能量E），非定常守恒型動(dòng)網(wǎng)格流場(chǎng)計(jì)算方程為：

式（6）中：V(t)為空間中大小和形狀都隨時(shí)間變化的控制體積；ρ為流體密度；?V( )t為控制體積的運(yùn)動(dòng)邊界；ug為運(yùn)動(dòng)網(wǎng)格的運(yùn)動(dòng)速度；u為流體速度矢量；Γ為耗散系數(shù)；Sφ是標(biāo)量φ的源項(xiàng)。

湍流模型采用計(jì)算精度比較高、應(yīng)用比較廣泛的k-ε二方程模型。

k控制方程：

ε控制方程：

式（7）（8）中：P為湍流動(dòng)能產(chǎn)生項(xiàng)；vt為黏性系數(shù)，vt=，而μt=；Cμ、Cε1、Cε2、σε和σk為模型系數(shù)。

另外，為使方程組封閉，還有氣體的狀態(tài)方程：

式（9）中：pg為氣體壓力；R為氣體常數(shù)；T為氣體溫度，單位K。

2 數(shù)值仿真方法

2.1 流場(chǎng)網(wǎng)格離散模型

本文為簡(jiǎn)化計(jì)算，僅考慮炸彈質(zhì)心運(yùn)動(dòng)位置及俯仰角度變化，不考慮炸彈滾動(dòng)、偏航角度，因此，這里采用二維流場(chǎng)仿真計(jì)算。選擇飛機(jī)艙彈分離的炸彈部分運(yùn)動(dòng)區(qū)域作為仿真區(qū)域，采用有限體積法[17]的非結(jié)構(gòu)化網(wǎng)格離散這個(gè)區(qū)域，最終所建立的流場(chǎng)仿真區(qū)域離散網(wǎng)格及邊界如圖2所示。將流場(chǎng)來(lái)流方向設(shè)置為壓力進(jìn)口邊界，出流方向設(shè)置為遠(yuǎn)場(chǎng)邊界。

圖2 流場(chǎng)離散網(wǎng)格Fig.2 Discrete grid of flow field

2.2 動(dòng)網(wǎng)格仿真技術(shù)

為保證流場(chǎng)仿真區(qū)域內(nèi)剛體（炸彈）氣動(dòng)力的計(jì)算精度，首先，在彈體周圍生成流體邊界層，讓流體邊界層隨著彈體一起運(yùn)動(dòng)，保證彈體周圍的邊界層不變，并且在彈體附近區(qū)域加密計(jì)算網(wǎng)格[15]。在計(jì)算過程中，為避免由于剛體運(yùn)動(dòng)導(dǎo)致流場(chǎng)區(qū)域網(wǎng)格扭曲，品質(zhì)變壞，嚴(yán)重影響仿真精度的情況發(fā)生，這里采用彈性光順法（smoothing）和局部網(wǎng)格重構(gòu)法（remeshing）2 項(xiàng)動(dòng)網(wǎng)格技術(shù)[17]。彈性光順法能夠保證整個(gè)仿真區(qū)域網(wǎng)格節(jié)點(diǎn)像彈簧連接的網(wǎng)格系統(tǒng)一樣，在計(jì)算時(shí)間步更新后，重新達(dá)到新的平衡位置，減少全域網(wǎng)格扭曲變形；局部網(wǎng)格重構(gòu)法保證了剛體附近局部網(wǎng)格扭曲率或尺寸超過設(shè)定標(biāo)準(zhǔn)時(shí)，局部網(wǎng)格將被重新劃分，從而減少了局部區(qū)域網(wǎng)格過大變形，但同時(shí)，局部網(wǎng)格間連接屬性、節(jié)點(diǎn)數(shù)量和連接關(guān)系發(fā)生改變。應(yīng)用彈性光順法和局部網(wǎng)格重構(gòu)法后，網(wǎng)格質(zhì)量得到顯著改善。流場(chǎng)仿真網(wǎng)格如圖3所示。

圖3 網(wǎng)格重構(gòu)后的流場(chǎng)仿真網(wǎng)格Fig.3 Flow field simulation network after grid reconstruction

2.3 仿真耦合方法

本文采用CFD和RBD方程耦合求解來(lái)仿真艙彈分離過程[18-19]：首先，利用CFD方程計(jì)算出某一時(shí)刻仿真區(qū)域流場(chǎng)參數(shù)分布，獲得流場(chǎng)中炸彈所受氣動(dòng)力和氣動(dòng)力矩，將其傳遞給RBD方程以獲得炸彈對(duì)氣動(dòng)力的響應(yīng)；然后，通過RBD 方程計(jì)算出下一時(shí)刻的炸彈運(yùn)動(dòng)位置和姿態(tài)，根據(jù)這些信息更新計(jì)算網(wǎng)格（經(jīng)過反復(fù)迭代直至滿足動(dòng)網(wǎng)格設(shè)置標(biāo)準(zhǔn)）；再進(jìn)行下一時(shí)刻的CFD方程計(jì)算以獲取新的氣動(dòng)力；重復(fù)以上耦合過程，直至計(jì)算完畢。

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

本文仿真的剛體對(duì)象——炸彈質(zhì)量m為500 kg，轉(zhuǎn)動(dòng)慣量Lzz為500 kg?m2，飛機(jī)飛行馬赫數(shù)為0.76，初始流場(chǎng)區(qū)域劃分的網(wǎng)格單元數(shù)為19 596，最小單元面積為1.827×10-2m2，最大單元面積為1.908 m2。在艙彈分離時(shí)刻，飛機(jī)向炸彈施加了拋射力（轉(zhuǎn)換到重心上為Fcx、Fcy、Mcz），作用時(shí)間0.3 s。在Fcx=-10 000 N、Fcy=-8 000 N 下，Mcz分別為0 N?m、2 200 N?m、-2 200 N?m，計(jì)算不同力矩下炸彈姿態(tài)變化。

當(dāng)Mcz=0 N?m 時(shí)，炸彈在拋射力Fc作用下加速離開彈艙向下運(yùn)動(dòng)，出艙瞬間在外界干擾下極易發(fā)生彈體反轉(zhuǎn)，如圖4所示。在機(jī)艙艙門剛打開，炸彈離開原位瞬間，由于彈體上下表面受外界氣流壓力不均衡，炸彈頭部微微向下擺動(dòng)（轉(zhuǎn)動(dòng)角度為2°），0.4 s 彈體完全出艙脫落飛機(jī)邊界層氣流的影響，在外界氣流作用下，彈體周圍壓力分布如圖5 所示。彈體上氣體壓力、黏性力合成出1 個(gè)順時(shí)針氣動(dòng)力矩推動(dòng)彈體開始做順時(shí)針轉(zhuǎn)動(dòng)，彈體轉(zhuǎn)動(dòng)角度由正值變?yōu)樨?fù)值，如圖6 所示。在1.38 s 時(shí)轉(zhuǎn)動(dòng)角度達(dá)到-90°，頭部朝上，轉(zhuǎn)動(dòng)速度仍為負(fù)值，彈體在轉(zhuǎn)動(dòng)慣性作用下姿態(tài)角度繼續(xù)增大。此后，炸彈在重力、氣動(dòng)力作用下，彈體姿態(tài)角度迅速調(diào)整到-270°，彈頭向下。

圖4 炸彈下落軌跡及姿態(tài)（Mcz=0 N?m）Fig.4 Falling track and attitude of bomb(Mcz=0 N?m)

圖5 彈體周圍壓力分布Fig.5 Pressure distribution around the bomb

圖6 彈體姿態(tài)角度變化Fig.6 Attitude angle change of bomb

當(dāng)炸彈出艙時(shí)向其施加1 個(gè)正向力矩作用（Mcz=2 200 N?m），彈體獲得了1 個(gè)逆時(shí)針轉(zhuǎn)動(dòng)角速度。炸彈出艙后，在重力、氣動(dòng)力矩的作用下，彈頭迅速向下擺動(dòng)，并且轉(zhuǎn)動(dòng)角速度也迅速增大，1.1 s 彈體轉(zhuǎn)動(dòng)角度就達(dá)到90°，彈頭擺正到垂直向下打擊姿態(tài)，如圖7所示。此后，彈體主要在慣性力矩作用下，圍繞垂直打擊姿態(tài)（炸彈軸線與重力作用線重合）做左右調(diào)整擺動(dòng)，打擊姿態(tài)迅速穩(wěn)定。

圖7 炸彈下落姿態(tài)（Mcz=2 200 N?m）Fig.7 Falling attitude of bomb(Mcz=2 200 N?m)

當(dāng)炸彈出艙時(shí)向其施加1 個(gè)負(fù)向力矩時(shí)（Mcz=-2 200 N?m），炸彈出艙后，彈體轉(zhuǎn)動(dòng)角度、轉(zhuǎn)動(dòng)角度順時(shí)針迅速增加，1.12 s 時(shí)轉(zhuǎn)動(dòng)角度達(dá)-90°，2.23 s時(shí)彈體轉(zhuǎn)動(dòng)到-270°，導(dǎo)彈垂直打擊姿態(tài)如圖8所示。

圖8 炸彈下落姿態(tài)（Mcz=-2 200 N?m）Fig.8 Falling attitude of bomb(Mcz=-2 200 N?m)

4 結(jié)論

本文基于網(wǎng)格彈性光順-重構(gòu)的動(dòng)網(wǎng)格技術(shù)，將CFD與RBD方程耦合，仿真計(jì)算了炸彈從飛機(jī)艙分離后姿態(tài)角度變化歷程，仿真過程及結(jié)果表明：

1）采用彈性光順法和局部網(wǎng)格重構(gòu)法相結(jié)合的動(dòng)網(wǎng)格技術(shù)，能夠有效地避免剛體運(yùn)動(dòng)引起的網(wǎng)格畸變，可顯著提高網(wǎng)格品質(zhì)和計(jì)算精度；

2）在艙彈分離時(shí)，向炸彈施加不同的力矩，彈體姿態(tài)調(diào)整方向是不一樣的，當(dāng)向其施加1 個(gè)正向的拋射力矩，有助于炸彈在下降過程中快速地?cái)[正到垂直打擊姿態(tài)，縮短姿態(tài)調(diào)整時(shí)間。