基于計算流體力學(xué)與剛體動力學(xué)耦合的高速旋轉(zhuǎn)彈丸彈道計算方法

鐘陽，王良明，吳映鋒,2

(1.南京理工大學(xué) 能源與動力工程學(xué)院, 江蘇 南京 210094; 2.63961 部隊, 北京 100012)

0 引言

先進(jìn)的武器裝備是國防建設(shè)的必要條件，射程、精度、威力等是武器系統(tǒng)的主要戰(zhàn)術(shù)技術(shù)指標(biāo)，武器系統(tǒng)的研制常從彈道設(shè)計[1]開始，因此彈道計算在整個武器系統(tǒng)研發(fā)過程中有著非常重要的地位。傳統(tǒng)的彈道計算方法是利用基于氣動力和力矩系數(shù)的彈道模型，通過龍格- 庫塔(Runge-Kutta)法或阿當(dāng)姆斯預(yù)報- 校正法進(jìn)行數(shù)值求解。傳統(tǒng)彈道計算方法需要建立彈箭氣動力模型，對于旋成體彈丸，氣動模型相對準(zhǔn)確；而對于氣動外形復(fù)雜甚至非對稱的彈箭，建立準(zhǔn)確的氣動模型非常困難。因此，對于外形復(fù)雜的彈丸，如何準(zhǔn)確獲取其動態(tài)變化過程和彈道飛行特性，是彈道學(xué)所面臨的一個難題。

近十多年發(fā)展起來的計算流體力學(xué)(CFD)和剛體動力學(xué)(RBD)耦合計算方法，為彈道仿真提供了一種新的技術(shù)途徑。相對于傳統(tǒng)方法，該方法能夠?qū)崟r模擬彈箭在大氣中真實(shí)的飛行情況，即使再復(fù)雜的彈箭，也能從實(shí)時變化的流場參數(shù)中獲取氣動力和力矩，在無需提供氣動模型情況下完成彈道計算。另外，該方法能夠?qū)崟r提供彈箭飛行過程中彈體表面和周圍空氣的壓力、密度、溫度等全息流場參數(shù)，為研究人員提供豐富的數(shù)據(jù)資料，有利于新型彈箭的研制。文獻(xiàn)[2]通過歐拉角建立彈體坐標(biāo)系與地面坐標(biāo)系的關(guān)系，在彈體坐標(biāo)系下采用任意拉格朗日- 歐拉(ALE)控制方程對流場進(jìn)行求解，首次將CFD/RBD耦合技術(shù)成功應(yīng)用到彈道計算中。進(jìn)而以某尾翼穩(wěn)定彈為對象，采用阿當(dāng)姆斯預(yù)報- 校正法對流場和6自由度剛體動力學(xué)方程組進(jìn)行耦合求解，并通過試驗(yàn)結(jié)果驗(yàn)證了方法的有效性。該方法被進(jìn)一步應(yīng)用于超音速、跨音速和亞音速彈丸模擬中，通過彈道仿真結(jié)果和試驗(yàn)結(jié)果對比再一次驗(yàn)證了方法的正確性[3-5]。文獻(xiàn)[6]在耦合CFD/RBD技術(shù)中引入入口邊界條件，研究了噴流對彈道的影響規(guī)律，驗(yàn)證了在耦合計算中加入進(jìn)口邊界來模擬脈沖修正彈彈道的可行性。文獻(xiàn)[7-11]通過彈丸飛行動力學(xué)模型推導(dǎo)出氣動參數(shù)估計模型，根據(jù)耦合CFD/RBD計算結(jié)果估計出彈丸的氣動系數(shù)。文獻(xiàn)[12-16]在耦合CFD/RBD技術(shù)中加入飛行控制系統(tǒng)(FCS)，建立耦合CFD/RBD/FCS方法并模擬了鴨舵受控尾翼彈的姿態(tài)運(yùn)動，通過與實(shí)驗(yàn)數(shù)據(jù)對比驗(yàn)證了方法的正確性。文獻(xiàn)[17]采用CFD/RBD和非結(jié)構(gòu)動網(wǎng)格技術(shù)計算均勻流場下不同密度平板的運(yùn)動軌跡，分析了碎片對飛機(jī)造成的二次傷害。文獻(xiàn)[18]基于自主開發(fā)軟件HUNS3D，結(jié)合非結(jié)構(gòu)嵌套網(wǎng)格技術(shù)，使用阿當(dāng)姆斯預(yù)估- 校正法耦合求解6自由度剛體運(yùn)動方程，研究了子母彈分離過程。文獻(xiàn)[19]采用剛性動網(wǎng)格和CFD技術(shù)，在規(guī)定運(yùn)動下模擬了帶鴨舵彈丸的運(yùn)動過程，研究了錐形運(yùn)動下彈丸非定常空氣動力學(xué)特性。

綜上所述，文獻(xiàn)[19]采用的是規(guī)定彈丸運(yùn)動的方法，不能直接用于研究彈丸的自由運(yùn)動。文獻(xiàn)[17-18]這類方法的計算域通常包含物體運(yùn)動軌跡，并不適用于研究射程達(dá)幾十千米的彈丸。相比較而言，文獻(xiàn)[2-16]基于彈體系建立流動關(guān)系更為合理。然而，對于高速旋轉(zhuǎn)彈丸，基于彈體系建立流動關(guān)系也存在不足之處，高速旋轉(zhuǎn)彈丸會使遠(yuǎn)場產(chǎn)生很大的附加速度，造成遠(yuǎn)場計算不準(zhǔn)確[20]?，F(xiàn)有的CFD/RBD耦合方法一般采用阿當(dāng)姆斯預(yù)報- 校正法，該方法較容易實(shí)現(xiàn)，但不能自啟動[1]。在流場初始耦合時只能采用當(dāng)前值來代替過去值，這必將在初始階段損失時間精度，且不易變步長。

本文針對已有研究的不足，基于ALE形式的流動模型，建立控制體表面運(yùn)動和彈軸坐標(biāo)系運(yùn)動的耦合模型。在此基礎(chǔ)上提出一種基于4階Runge-Kutta法的CFD/RBD耦合計算方法，為研究高速旋轉(zhuǎn)彈丸的真實(shí)彈道和流場提供參考。

1 CFD數(shù)值計算方法

1.1 基于彈軸運(yùn)動的ALE形式控制方程

彈丸飛行時姿態(tài)不斷變化，與彈丸固連的計算網(wǎng)格具有附加速度，ALE形式的Navier-Stokes方程(簡稱N-S方程)能夠有效地求解具有網(wǎng)格速度的流場。由于彈丸高速旋轉(zhuǎn)，在彈體系下建立流動關(guān)系時，遠(yuǎn)場流動速度相對于壁面附近非常大，從而導(dǎo)致計算不準(zhǔn)確[20]。而彈軸坐標(biāo)系有效去除了彈體高速自轉(zhuǎn)，克服了這一問題?；趶椵S運(yùn)動建立具有ALE形式的N-S方程并對流動進(jìn)行求解，其積分形式為

(1)

(2)

Fc為靜態(tài)網(wǎng)格下的對流通量，v?Ω為附加在控制體表面的反變速度。

當(dāng)彈軸系的平動速度矢量和轉(zhuǎn)動角速度矢量分別為vp和ωa時，附加在控制體表面的反變速度為

v?Ω=n{vp+ωa×(r-rg)}，

(3)

式中：n=[nx,ny,nz]T為控制體表面的單位外法向矢量，nx、ny、nz分別為單位外法向矢量在x軸、y軸、z軸上的分量；r為控制體表面微元的矢徑；rg為彈丸質(zhì)心矢徑。

1.2 修正簡單低耗散迎風(fēng)矢通量分裂格式

(1)式中對流通量的計算方法是CFD的研究熱點(diǎn)之一。簡單低耗散迎風(fēng)矢通量分裂格式(SLAU)是著名迎風(fēng)型矢通量分裂格式(AUSM)系列的一個發(fā)展方向，具有低耗散、格式簡潔、低馬赫數(shù)下無需調(diào)整參數(shù)等優(yōu)點(diǎn)，其改進(jìn)型SLAU2[21]在壓力通量中加入了與來流馬赫數(shù)相關(guān)的耗散項(xiàng)，繼承了原格式優(yōu)點(diǎn)，同時實(shí)現(xiàn)了比其他全速格式更簡潔的形式。引入彈丸運(yùn)動后，基于SLAU2的計算格式為

(4)

(5)

cij=(cl+cr)/2,

(6)

式中：cl和cr分別為控制面左右兩側(cè)音速。

壓力通量為

(7)

(8)

(9)

M為法向流動速度與音速比值，M=(v·n)/cij.

1.3 邊界條件、時間推進(jìn)及湍流模型

(10)

(11)

式中：下標(biāo)i表示控制體第i個面；nf為面的個數(shù)；vi為控制體第i個面上的流動速度矢量；ni為控制體第i個面的單位外法向矢量；ci和ΔSi分別為離散后控制體第i個面的音速和面積。

計算(1)式中黏性通量的湍流黏性時采用Spalart-Allmaras(S-A)湍流模型，該模型具有計算量小、穩(wěn)定性好、計算精度滿足工程要求等優(yōu)點(diǎn)，在氣動力計算和航空航天等領(lǐng)域得到廣泛應(yīng)用。S-A模型發(fā)展了較多版本，本文通過某彈丸三維外部繞流數(shù)值實(shí)驗(yàn)發(fā)現(xiàn)，標(biāo)準(zhǔn)可壓S-A湍流模型[23]對于高雷諾數(shù)和低雷諾數(shù)的適應(yīng)能力較強(qiáng)，且在尾部網(wǎng)格不對稱情況下能夠收斂到較好的結(jié)果，因此本文采用該湍流模型計算湍流黏性系數(shù)。

2 基于CFD的彈丸飛行動力學(xué)模型

CFD坐標(biāo)系通常以彈頭為原點(diǎn)，Oxc軸與彈軸重合指向彈尾為正，Oyc軸垂直O(jiān)xc軸向上為正，Ozc軸滿足右手法則。其他坐標(biāo)系無特殊說明均參考文獻(xiàn)[1]，彈軸坐標(biāo)系無特殊說明均指第1彈軸坐標(biāo)系。旋轉(zhuǎn)穩(wěn)定彈丸質(zhì)心運(yùn)動的動力學(xué)方程組通常采用彈道坐標(biāo)系Otxtytzt下的力系形式[1]。為了使CFD計算的空氣動力和空氣動力矩能夠耦合到彈道模型中，需建立轉(zhuǎn)換關(guān)系。首先將CFD坐標(biāo)系下的力系轉(zhuǎn)換到第1彈軸坐標(biāo)系下，再將第1彈軸坐標(biāo)系下力系轉(zhuǎn)到第2彈軸坐標(biāo)系下，最后將第2彈軸坐標(biāo)系下的力系轉(zhuǎn)到彈道坐標(biāo)系下。CFD坐標(biāo)系到第1彈軸坐標(biāo)系的轉(zhuǎn)換矩陣為

Lac=diag[-1 1 -1].

(12)

第1彈軸坐標(biāo)系到彈道坐標(biāo)系的轉(zhuǎn)換矩陣為

(13)

(14)

式中：vt為彈丸質(zhì)心運(yùn)動速度；θa為速度高低角；ψd為速度方向角；φd為彈軸方位角；m為彈丸質(zhì)量；g為重力加速度；Fcξ、Fcη、Fcζ分別為空氣動力在CFD坐標(biāo)系Oxc軸、Oyc軸、Ozc軸上的分量；Lc為系數(shù)矩陣，

Lc=diag[1/m1/mvtcosψd1/mvt].

(15)

彈丸質(zhì)心運(yùn)動的運(yùn)動學(xué)方程組[1]為

(16)

式中：vpx、vpy、vpz分別為彈丸質(zhì)心速度矢量vp在地面坐標(biāo)系Oexeyeze中Oexe軸、Oeye軸、Oeze軸上的分量。由(12)式和基準(zhǔn)坐標(biāo)系到彈軸坐標(biāo)系的轉(zhuǎn)換矩陣[1]，得到控制體表面附加平動速度分量為

(17)

式中：vcξ、vcη、vcζ分別為控制體表面附加平動速度矢量；φa為彈軸高低角。

根據(jù)(12)式，得到彈丸總空氣動力矩矢量M在彈軸坐標(biāo)系下的3個分量為

(18)

式中：Mcξ、Mcη、Mcζ分別為CFD數(shù)值計算獲得的彈丸總空氣動力矩矢量在CFD坐標(biāo)系3軸上的投影。將(18)式代入彈丸在彈軸坐標(biāo)系下繞質(zhì)心轉(zhuǎn)動的動力學(xué)方程組[1]，可得

(19)

式中：C為極轉(zhuǎn)動慣量；A為赤道轉(zhuǎn)動慣量；ωξ、ωη和ωζ分別為彈丸角速度矢量在彈軸坐標(biāo)系上3個軸上的投影分量。由(12)式得到壁面角速度為

(ωcwξ,ωcwη,ωcwζ)=(-ωξ,ωη,-ωζ).

(20)

彈軸坐標(biāo)系角速度矢量ωa在彈軸坐標(biāo)系3軸分量為(ωaξ,ωaη,ωaζ)=(ωζtanφd,ωη,ωζ)[1]，通過轉(zhuǎn)換矩陣(12)式得到控制體表面附加角速度為

(ωcξ,ωcη,ωcζ)=(-ωζtanφd,ωη,-ωζ)，

(21)

對于高速旋轉(zhuǎn)彈丸，由于|ωζtanφd|遠(yuǎn)小于|ωξ|，由(3)式可知，在彈軸坐標(biāo)系下控制體表面附加速度遠(yuǎn)小于彈體坐標(biāo)系，特別是遠(yuǎn)場。

彈道方程組中繞質(zhì)心轉(zhuǎn)動的運(yùn)動學(xué)方程組和δa、δd、β3個角關(guān)系方程參見文獻(xiàn)[1]，不再贅述。

3 CFD/RBD耦合計算方法

本文彈道和流場耦合計算框架如圖1所示，本節(jié)主要研究圖1中的耦合模塊。

圖1 CFD/RBD耦合計算框架Fig.1 Framework of CFD/RBD coupling calculation

(14)式、(16)式、(19)式和文獻(xiàn)[1]中的繞心轉(zhuǎn)動運(yùn)動學(xué)方程組以及3個角關(guān)系方程，共同組成了由CFD數(shù)值計算的空氣動力和空氣動力矩驅(qū)動彈道方程。(17)式、(20)式和(21)式將彈道方程解算出的狀態(tài)量反饋給CFD求解器，從而形成了耦合機(jī)制。為敘述方便，(17)式、(20)式和(21)式統(tǒng)稱為接口模型。

耦合計算采用經(jīng)典的4階Runge-Kutta法。已知n時間層面值(tn,y1,n,y2,n,…,ym,n)，tn表示數(shù)值計算到第n步的時間，ym,n表示第n步下第m個參數(shù)，步長為h. 則求解n+1時間層面各函數(shù)值的公式為

(22)

式中：

(23)

fi表示第i個微分方程的右端表達(dá)式。由(23)式可知，每一個Runge-Kutta子步都需要?dú)鈩訑?shù)據(jù)的支持。有2種較簡單的氣動力和力矩耦合方法：第1種為定值法，借鑒系數(shù)凍結(jié)法的思想，將n時間層面CFD計算出的瞬態(tài)氣動力和力矩凍結(jié)，即在Runge-Kutta 4個子步中保持不變；第2種為插值法，利用n-1和n時間層面的氣動力和力矩，在4個子步中進(jìn)行線性外插。本文提出一種緊耦合實(shí)現(xiàn)方法，使每一個子步中的氣動力和力矩都為CFD計算出的瞬態(tài)結(jié)果。具體實(shí)現(xiàn)步驟如下：

步驟1保存n時間層面的流場參數(shù)，并積分出CFD坐標(biāo)系下彈丸受到的空氣動力和力矩，通過彈道方程右端子式計算出ki1.

步驟2將流場參數(shù)設(shè)置為步驟1保存的n時間層面；根據(jù)ki2右端括號中的參數(shù)形式，利用步驟1計算出的ki1和接口模型計算出控制體表面附加運(yùn)動速度。將時間步長設(shè)為0.5h，采用雙時間步將流場推進(jìn)到tn+0.5h時刻；積分出彈丸受到的空氣動力和力矩。利用彈道方程右端表達(dá)式計算出ki2.

步驟3將流動情況返回到n時間層面；根據(jù)ki3右端參數(shù)形式，利用步驟2計算出的ki2和接口模型計算出控制體表面附加運(yùn)動速度；雙時間步步長設(shè)置為0.5h將流場推進(jìn)到tn+0.5h時刻；積分出作用在彈丸上的空氣動力和力矩；利用彈道方程右端表達(dá)式計算出ki3.

步驟4將流動數(shù)設(shè)置為tn時刻的值，通過步驟3計算出的ki3來獲得ki4右端參數(shù)值，由接口模型獲取控制體表面附加運(yùn)動情況；以h為步長，采用雙時間步將流場推進(jìn)到tn+h時刻并獲取彈丸上受到的氣動力和力矩；由彈道方程右端表達(dá)式計算出ki4.

步驟5由(22)式計算出n+1時間層面的彈道參數(shù)，并通過接口模型計算出控制體表面附加運(yùn)動；將流場返回至tn時刻，以h為步長，利用雙時間步將流場推進(jìn)至tn+h時刻；積分出當(dāng)前的空氣動力和力矩。

圖2 流場和彈道緊耦合計算過程Fig.2 Flow field and ballistic tight coupling calculation process

通過以上5步可完成流場和彈道方程從n時間層面到n+1時間層面完整的4階Runge-Kutta緊耦合計算過程，計算流程示意圖如圖2所示，圖2中j=1，2，3，4表示Runge-Kutta4個步驟。需要注意的是，步驟2～步驟4中涉及到的tn+0.5h和tn+h時刻只是中間過程，步驟5中采用雙時間步推進(jìn)到的tn+h時刻為n+1時間層面。

4 仿真驗(yàn)證與分析

4.1 M549彈丸定態(tài)運(yùn)動和規(guī)定運(yùn)動模擬

計算模型采用美國M549旋成體彈丸，具體結(jié)構(gòu)尺寸參見文獻(xiàn)[24]。利用Solidworks商業(yè)軟件進(jìn)行幾何建模，如圖3所示。

圖3 M549彈丸模型Fig.3 M549 projectile model

將幾何模型導(dǎo)入CFD前處理軟件ICEM中進(jìn)行網(wǎng)格劃分，如圖4所示。文獻(xiàn)[25]對旋成體彈丸的網(wǎng)格劃分進(jìn)行了詳細(xì)研究，并通過與試驗(yàn)數(shù)據(jù)對比驗(yàn)證了網(wǎng)格劃分的合理性。本文在網(wǎng)格劃分過程中部分技術(shù)參數(shù)參考此文獻(xiàn)，如遠(yuǎn)場位置、近壁第1層網(wǎng)格厚度和黏性層網(wǎng)格延展率等。為使彈丸頭部和尾部網(wǎng)格具有良好的正交性，進(jìn)行O型剖分，網(wǎng)格總數(shù)約為186萬。

圖4 M549彈丸計算域網(wǎng)格Fig.4 Computational domain mesh of M549 projectile

通過使彈丸做定態(tài)飛行和規(guī)定運(yùn)動獲得彈丸的氣動特性，將結(jié)果與Spinner數(shù)據(jù)[24]進(jìn)行對比，以檢驗(yàn)本文CFD數(shù)值方法的可靠性。

(24)

圖5 俯仰力矩系數(shù)遲滯回線Fig.5 Hysteresis loop of pitching moment coefficient

表1 馬赫數(shù)為2時M549彈丸氣動系數(shù)Tab.1 Aerodynamic coefficients of M549 projectile for Ma=2

總之，在馬赫數(shù)為2.0時，靜態(tài)氣動系數(shù)和動態(tài)氣動系數(shù)的仿真結(jié)果相對Spinner數(shù)據(jù)的誤差分別小于6%和10%，數(shù)值仿真的氣動系數(shù)較可靠。

4.2 M549彈丸自由飛行仿真驗(yàn)證與分析

文獻(xiàn)[27]基于傳統(tǒng)彈道計算方法，采用旋成體彈丸氣動模型和CFD仿真氣動系數(shù)計算某旋成體彈丸彈道，并與靶道試驗(yàn)結(jié)果對比，驗(yàn)證了CFD數(shù)值方法結(jié)合傳統(tǒng)彈道計算方法預(yù)測旋成體彈丸實(shí)際彈道的可行性。本文借鑒文獻(xiàn)[27]的方法計算M549旋成體彈丸彈道，并以該彈道結(jié)果為基準(zhǔn)，研究第3節(jié)中3種力和力矩耦合方法以及不同時間步長對耦合計算方法模擬彈道能力的影響。

彈道起始條件為：射角45°；射向0°；高低攻角和飛行馬赫數(shù)分別沿用4.1節(jié)中的3°和2.0；方向攻角0°；彈軸無擺動。流場仿真條件為：雙時間步子迭代最大計算步數(shù)為200；收斂條件為10-2；流場初始化采用4.1節(jié)中對應(yīng)條件下的計算收斂結(jié)果。

對于傳統(tǒng)6自由度彈道計算，時間步長通常取5 ms左右[1]。為了研究3種力和力矩耦合方法求解彈道的能力，取步長h=5 ms進(jìn)行仿真，得到攻角關(guān)系曲線如圖6所示。通過觀察圖6中的攻角變化情況發(fā)現(xiàn)：采用定值法效果最差，仿真攻角結(jié)果發(fā)散，計算彈道失??；插值法和緊耦合法均模擬出一般旋轉(zhuǎn)穩(wěn)定彈丸的雙圓運(yùn)動，其中緊耦合法計算結(jié)果與基準(zhǔn)結(jié)果更加接近，即緊耦合精度最好。

圖6 h=5 ms時3種耦合方法對比Fig.6 Comparison of 3 coupling methods for h=5 ms

對于非定常流場，為了準(zhǔn)確反映流動細(xì)節(jié)，時間步長要小得多。流場與彈道耦合計算的主要矛盾是非定常氣動參數(shù)的收斂程度，對某155 mm口徑彈丸進(jìn)行非定常氣動力計算時發(fā)現(xiàn)，時間步長取0.5 ms和0.05 ms時滾轉(zhuǎn)力矩系數(shù)差別僅為4.67%[28]，取時間步長為0.5 ms時兼具精度和效率。

圖7 彈道參數(shù)仿真結(jié)果Fig.7 Simulated results of trajectory parameters

采用緊耦合方法，在時間步長分別為5 ms和0.5 ms條件下進(jìn)行仿真，仿真結(jié)果與基準(zhǔn)結(jié)果如圖7所示。由圖7(a)可見：總攻角δt初始為3°，在陀螺穩(wěn)定力矩作用下逐漸衰減；步長5 ms下總攻角計算結(jié)果偏大，衰減明顯不足；步長0.5 ms下衰減明顯加快，且衰減速度能夠跟隨基準(zhǔn)結(jié)果；從二者均值上看，初值均在4°左右，當(dāng)時間為1.805 s時，前者衰減到1.714 8°，后者衰減到1.170 4°，小步長比大步長衰減效率提高了約24%. 從圖7(b)和圖7(c)中彈軸擺動曲線以及圖7(e)中速度方向角曲線可以看出：當(dāng)步長為5 ms時仿真結(jié)果在第1個快圓運(yùn)動周期內(nèi)與基準(zhǔn)結(jié)果基本一致，但之后不能保持，雖然曲線形狀基本一致，但幅值偏大；當(dāng)步長為0.5 ms時，仿真結(jié)果與基準(zhǔn)結(jié)果始終保持較高的符合度。小步長下速度高低角與基準(zhǔn)結(jié)果基本重合，大步長下在波谷處略低，如圖7(d)所示。彈道側(cè)偏ze反映了地面坐標(biāo)系Oexeyeze中彈丸偏離射擊面Oexeye[1]的程度，如圖7(f)所示，仿真結(jié)果與基準(zhǔn)結(jié)果符合較好，在t=1.80 s時，0.5 ms步長和5 ms步長下的彈道側(cè)偏量與基準(zhǔn)結(jié)果偏差分別為2.18%和4.92%.

4.3 修正組件不轉(zhuǎn)時二維彈道修正彈模擬與分析

旋轉(zhuǎn)穩(wěn)定二維彈道修正彈通過后體高速旋轉(zhuǎn)來保持飛行穩(wěn)定?？刂茣r，修正組件滾轉(zhuǎn)角固定在某個位置，對彈道進(jìn)行修正。此時亦可采用本文方法進(jìn)行仿真。計算模型和網(wǎng)格劃分如圖8所示，網(wǎng)格劃分技術(shù)參數(shù)與4.1節(jié)相同，網(wǎng)格數(shù)約260萬。仿真時，修正組件滾轉(zhuǎn)角固定在圖8(a)位置即滾轉(zhuǎn)角為0°；初始高低攻角為0°；后體初始轉(zhuǎn)速1 500 rad/s；其余條件與4.2節(jié)相同。圖8中：F1和F3表示一對差動舵受到的空氣動力；F2和F4表示一對控制舵受到的空氣動力。

圖8 二維彈道修正彈計算模型和網(wǎng)格Fig.8 Calculation model and meshes of two-dimensional trajectory correction projectile

圖9給出了步長0.5 ms下緊耦合法計算的攻角曲線及特征點(diǎn)P1～P6，攻角變化規(guī)律為雙圓運(yùn)動，并與解析方法[29]計算的結(jié)果對比，兩種方法在頻率和幅值方面大致相似，但不完全相同。從攻角收斂趨勢來看，控制力向上舵偏角產(chǎn)生的攻角向下，與鴨舵布局的尾翼穩(wěn)定彈相反，這一結(jié)論與文獻(xiàn)[29]一致。

圖9 攻角仿真結(jié)果Fig.9 Simulated result of angle of attack

圖10 彈丸表面壓力分布Fig.10 Pressure distribution on the surface of projectile

為了分析彈丸飛行過程中表面壓力變化情況，觀察特征點(diǎn)處的流場，如圖10所示，視圖方向垂直于彈道坐標(biāo)系中的Otxtyt平面，特征點(diǎn)參數(shù)如表2所示。由起始階段圖10(a)可以看出，與超音速來流直接接觸的頭部和鴨舵稍部壓力較高，最大處超過0.17 MPa；前后體連接處壓力分布在0.112～0.156 MPa之間，其中鴨舵正后方處在鴨舵激波內(nèi)部，壓力比周圍高；圖10(b)中彈軸擺動到了左上方，彈頭部壓力出現(xiàn)不對稱，迎風(fēng)區(qū)壓力增大；隨著彈軸從左上方擺到左下方，彈頭迎風(fēng)區(qū)變成右上方，且由于低頭使得水平舵有效迎角減小，壓力降低，如圖10(c)所示；當(dāng)彈軸從左下方位置擺到圖10(d)中下方略偏右位置時，上方差動舵有效迎角變大同時處于迎風(fēng)區(qū)，下方有效迎角減小同時處于背風(fēng)區(qū)，使得前者壓力增大后者壓力減?。划?dāng)彈軸從圖10(d)位置略微上擺至圖10(e)位置時，水平舵有效迎角增大，壓力也增大，直至擺動到圖10(f)時，彈丸近似完成一個慢圓周期。緊耦合計算中，隨著彈軸擺動，彈丸表面迎風(fēng)區(qū)域和背風(fēng)區(qū)域不斷變化，形成時變且復(fù)雜的壓力分布，壓力分布進(jìn)而又影響彈軸擺動，這與彈丸實(shí)際飛行過程是類似的。綜上所述可知，與普通彈道計算方法相比，緊耦合法計算不僅能向設(shè)計人員提供彈道數(shù)據(jù)，還能提供流場參數(shù)，展現(xiàn)更加實(shí)際的彈丸飛行過程。

表2 二維彈道修正彈攻角仿真結(jié)果Tab.2 Simulated results of angle of attack of two-dimensional trajectory correction projectile

5 結(jié)論

耦合計算流體力學(xué)和剛體動力學(xué)方法是外彈道學(xué)一個新的發(fā)展方向。本文考慮高速旋轉(zhuǎn)彈丸的特點(diǎn)，利用壁面速度表征旋成體表面的高速旋轉(zhuǎn)，結(jié)合高速旋轉(zhuǎn)彈丸彈道模型，將彈軸坐標(biāo)系平動速度和轉(zhuǎn)動角速度附加于控制體表面，建立了基于ALE形式新的流動模型。在此基礎(chǔ)上提出一種適合于高速旋轉(zhuǎn)彈丸的基于4階Runge-Kutta法流場和彈道耦合數(shù)值計算方法。最后通過對比傳統(tǒng)方法的仿真結(jié)果和解析結(jié)果驗(yàn)證了緊耦合計算方法的有效性。得到主要結(jié)論如下：

1)在進(jìn)行耦合計算時，氣動力和力矩耦合模型對計算結(jié)果影響較大，采用線性插值和提供實(shí)時瞬態(tài)值即緊耦合時計算結(jié)果收斂性好，精度高。

2)時間步長的選取很重要，在合理時間步長下進(jìn)行緊耦合計算時，可以有效模擬彈丸飛行過程。

3)緊耦合方法模擬具有復(fù)雜外形的旋轉(zhuǎn)穩(wěn)定彈飛行時，能夠反映出彈丸的運(yùn)動特性，流場變化規(guī)律符合客觀實(shí)際。

本文提出的耦合計算方法可為新型旋轉(zhuǎn)穩(wěn)定彈，如加裝了精確制導(dǎo)組件的二維彈道修正彈設(shè)計和計算提供一種新的理論和方法。今后可進(jìn)一步考慮控制部件的運(yùn)動，為新型旋轉(zhuǎn)穩(wěn)定彈控制組件的設(shè)計提供基礎(chǔ)。