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

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

0 引言

1 仿真計算方程

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

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

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

1.1 RBD方程

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

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

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

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

體坐標(biāo)系下炸彈角運動方程為：

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

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

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

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

1.2 CFD方程

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

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

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

k控制方程：

ε控制方程：

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

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

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

2 數(shù)值仿真方法

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

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

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

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

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

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

2.3 仿真耦合方法

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

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

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

當(dāng)Mcz=0 N?m 時，炸彈在拋射力Fc作用下加速離開彈艙向下運動，出艙瞬間在外界干擾下極易發(fā)生彈體反轉(zhuǎn)，如圖4所示。在機(jī)艙艙門剛打開，炸彈離開原位瞬間，由于彈體上下表面受外界氣流壓力不均衡，炸彈頭部微微向下擺動（轉(zhuǎn)動角度為2°），0.4 s 彈體完全出艙脫落飛機(jī)邊界層氣流的影響，在外界氣流作用下，彈體周圍壓力分布如圖5 所示。彈體上氣體壓力、黏性力合成出1 個順時針氣動力矩推動彈體開始做順時針轉(zhuǎn)動，彈體轉(zhuǎn)動角度由正值變?yōu)樨?fù)值，如圖6 所示。在1.38 s 時轉(zhuǎn)動角度達(dá)到-90°，頭部朝上，轉(zhuǎn)動速度仍為負(fù)值，彈體在轉(zhuǎn)動慣性作用下姿態(tài)角度繼續(xù)增大。此后，炸彈在重力、氣動力作用下，彈體姿態(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)炸彈出艙時向其施加1 個正向力矩作用（Mcz=2 200 N?m），彈體獲得了1 個逆時針轉(zhuǎn)動角速度。炸彈出艙后，在重力、氣動力矩的作用下，彈頭迅速向下擺動，并且轉(zhuǎn)動角速度也迅速增大，1.1 s 彈體轉(zhuǎn)動角度就達(dá)到90°，彈頭擺正到垂直向下打擊姿態(tài)，如圖7所示。此后，彈體主要在慣性力矩作用下，圍繞垂直打擊姿態(tài)（炸彈軸線與重力作用線重合）做左右調(diào)整擺動，打擊姿態(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)炸彈出艙時向其施加1 個負(fù)向力矩時（Mcz=-2 200 N?m），炸彈出艙后，彈體轉(zhuǎn)動角度、轉(zhuǎn)動角度順時針迅速增加，1.12 s 時轉(zhuǎn)動角度達(dá)-90°，2.23 s時彈體轉(zhuǎn)動到-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)的動網(wǎng)格技術(shù)，將CFD與RBD方程耦合，仿真計算了炸彈從飛機(jī)艙分離后姿態(tài)角度變化歷程，仿真過程及結(jié)果表明：

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

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