岳 通，王惠源，張成卿

(中北大學(xué) 機(jī)電工程學(xué)院，太原 030051)

彈箭氣動(dòng)參數(shù)對(duì)氣動(dòng)外形設(shè)計(jì)及彈道的計(jì)算具有重要意義，為此眾多學(xué)者對(duì)彈箭氣動(dòng)參數(shù)的計(jì)算進(jìn)行了大量的研究。文獻(xiàn)[1]研究了大長(zhǎng)徑比超音速卷弧尾翼火箭彈在標(biāo)準(zhǔn)氣象條件下,卷弧尾翼對(duì)全彈氣動(dòng)特性的影響；文獻(xiàn)[2]通過研究得到某35 mm超音速榴彈在標(biāo)準(zhǔn)氣象條件下阻力系數(shù)隨彈形的變化規(guī)律；文獻(xiàn)[3]研究了在標(biāo)準(zhǔn)氣象條件下，某型導(dǎo)彈的亞音速氣動(dòng)特性并給出了該彈的氣動(dòng)參數(shù)；文獻(xiàn)[4]研究了某小型單兵巡飛彈在2.5 km海拔高度處的氣動(dòng)參數(shù)，為巡飛彈的氣動(dòng)外形的設(shè)計(jì)提供參考；文獻(xiàn)[5]研究了 4 km高空條件下，有無彈帶對(duì)高速旋轉(zhuǎn)彈丸氣動(dòng)特性的影響，得到彈帶使得阻力系數(shù)有一定的升高，但彈帶形狀對(duì)彈丸升力系數(shù)影響不大的相關(guān)結(jié)論；文獻(xiàn)[6]采用計(jì)算流體動(dòng)力學(xué)對(duì)M910彈丸在全音速范圍內(nèi)靜壓為101 325 Pa,靜溫為 288 K條件下的氣動(dòng)特性進(jìn)行了模擬，得到各氣動(dòng)參數(shù)與實(shí)驗(yàn)數(shù)據(jù)的誤差在10%以內(nèi)，具有較高的計(jì)算精度。文獻(xiàn)[7]分析得到了某旋轉(zhuǎn)彈丸在全音速范圍內(nèi)，靜壓為101 325 Pa,靜溫為292 K條件下，靜態(tài)空氣動(dòng)力系數(shù)與大多數(shù)動(dòng)力導(dǎo)數(shù)具有很好的一致性。文獻(xiàn)[8]對(duì)來流馬赫數(shù)Ma=0.7～2.7，攻角為 0°情況下的旋轉(zhuǎn)彈丸繞流場(chǎng)進(jìn)行了數(shù)值模擬,所得阻力系數(shù)結(jié)果與實(shí)驗(yàn)值和半經(jīng)驗(yàn)公式計(jì)算值符合較好。文獻(xiàn)[9]通過運(yùn)用牛頓撞擊理論計(jì)算得到某高超聲速制導(dǎo)炮彈在靜壓為101 325 Pa，靜溫為293 K條件下的阻力系數(shù)與實(shí)驗(yàn)結(jié)果吻合較好。因此,該計(jì)算方法能夠可靠地獲得高超聲速制導(dǎo)炮彈的氣動(dòng)力特性。文獻(xiàn)[10]采用CFD方法，基于SST湍流模型對(duì)超聲速火箭彈在特定海拔高度下的外流場(chǎng)進(jìn)行了計(jì)算，計(jì)算結(jié)果表明該數(shù)值方法能較好地計(jì)算氣動(dòng)力參數(shù)。

然而在以往對(duì)彈箭氣動(dòng)參數(shù)的研究中，國內(nèi)外很多學(xué)者只研究了固定氣象條件下，彈箭的氣動(dòng)參數(shù)。少有人研究彈箭氣動(dòng)參數(shù)隨不同海拔高度處的變化規(guī)律，但是由于不同海拔高度處的氣象參數(shù)變化較大，則在不同海拔處發(fā)射的彈箭在飛行過程中受到的氣動(dòng)力也將會(huì)有較大的變化，從而影響彈箭的彈道。因此本文以小口徑尾翼彈為研究對(duì)象，研究該尾翼彈隨不同海拔高度，不同馬赫數(shù)處的氣動(dòng)特性及其外彈道特性的變化規(guī)律，為小口徑尾翼彈在不同海拔高度處的應(yīng)用奠定基礎(chǔ)。

1 數(shù)值計(jì)算方法

本文研究海拔高度在20 km以內(nèi)，小口徑尾翼彈的氣動(dòng)特性，其所屬流域?yàn)檫B續(xù)流域，氣體的間斷效應(yīng)并不顯著，連續(xù)性假設(shè)成立，所以仍可運(yùn)Navier-Stokes方程進(jìn)行計(jì)算，選用雷諾平均的N-S方程作為控制方程。

控制方程為[11]：

式中：U為解向量；F和G為通量向量；Ma∞為來流馬赫數(shù)；γ為比熱比，Re∞為來流雷諾數(shù)。

湍流模型采用渦黏模型中的Shear-Stress Transport(SST)k-ω模型，求解2個(gè)方程的輸運(yùn)方程，即湍動(dòng)能k,湍流頻率ω，SSTk-ω模型考慮到了湍流剪切應(yīng)力的輸運(yùn)，能較好的模擬彈箭的氣動(dòng)力，因此選用該湍流模型計(jì)算尾翼彈的氣動(dòng)參數(shù)，k和ω的輸運(yùn)方程具體形式為[12]：

尾翼彈外形參數(shù)如圖1所示。

圖1 尾翼彈外形參數(shù)

綜合考慮計(jì)算精度與計(jì)算效率，將計(jì)算域分層，在大的圓柱外流域內(nèi)再增加一個(gè)小圓柱加密區(qū)域[13]，大圓柱形外流域直徑為18倍彈徑長(zhǎng)度為10倍彈長(zhǎng)，小圓柱流域直徑為5倍彈徑，長(zhǎng)度為3倍彈長(zhǎng)。采用求解黏性子層的方法對(duì)近壁面區(qū)域進(jìn)行計(jì)算，即需要保證y+<5使得第一層網(wǎng)格節(jié)點(diǎn)位于粘性子層內(nèi)[14]，并保證邊界層數(shù)為8，劃分得到如圖2和圖3所示的多面體網(wǎng)格。

圖2 計(jì)算域網(wǎng)格

圖3 彈體表面網(wǎng)格

由美國1976年標(biāo)準(zhǔn)大氣分段逼近公式[15]計(jì)算：

計(jì)算位勢(shì)高度H,它與幾何高度Z的關(guān)系為：

H=Z/(1+Z/R0)

式中：R0=6356.766 km,為地球平均半徑；Z單位為 km。

1) 當(dāng)0≤Z≤11.019 km時(shí)：

p/p0=W5.252 9,ρ/ρ0=W4.255 9

2) 當(dāng)11.0191 km≤Z≤20.063 1 km時(shí)：

p/p0=0.119 53W，ρ/ρ0=0.158 98W

式中：p0=1 013.25 hPa，密度ρ0=1.225 kg/m3。

3) 黏性系數(shù)隨高度變化的計(jì)算公式：

式中：βa=1.458×10-6kg/(s·m·K1/2)，Ts=110.4 K。

式中：k為絕熱指數(shù)，空氣k=1.4;Rg為氣體常數(shù)，空氣Rg=287 J/(kg·K)；T為熱力學(xué)溫度(K)。

由以上公式計(jì)算得到各海拔高度的氣象參數(shù)值，作出相關(guān)曲線如圖4所示。

圖4 氣象參數(shù)變化曲線

設(shè)置流場(chǎng)外壁面為壓力遠(yuǎn)場(chǎng)邊界，彈體表面采用無滑移壁面邊界條件，求解方法采用壓力-速度耦合(Coupled)算法，對(duì)流通量采用二階迎風(fēng)格式進(jìn)行離散。對(duì)氣動(dòng)參數(shù)計(jì)算時(shí)，遠(yuǎn)場(chǎng)的溫度、壓力、密度、黏性系數(shù)等氣象參數(shù)已由以上部分計(jì)算得到。假設(shè)來流為理想氣體，來流馬赫數(shù)為Ma=0.6～3，攻角為α=0°。

算法有效性驗(yàn)證：

本文選用標(biāo)準(zhǔn)模型(Non-Rolling Basic Finner)進(jìn)行數(shù)值算法有效性驗(yàn)證，該標(biāo)準(zhǔn)模型已由美軍彈道研究實(shí)驗(yàn)室進(jìn)行了大量風(fēng)洞試驗(yàn)。圖5為其幾何外形參數(shù)，彈體直徑d=20 mm，以彈體橫截面積為參考面積,并用以上數(shù)值方法模擬標(biāo)模彈的外流場(chǎng)，得到零升阻力系數(shù)，然后作出圖6。

圖5 標(biāo)準(zhǔn)模型

圖6 標(biāo)模彈零升阻力系數(shù)曲線

圖6所示的仿真零升阻力系數(shù)曲線和文獻(xiàn)[16]的實(shí)驗(yàn)值擬合曲線大致吻合，誤差小于10%，在允許范圍以內(nèi)。由此可知，本文中所用的數(shù)值計(jì)算方法具有較高的可信度。

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

由數(shù)值方法計(jì)算得到不同海拔高度處的尾翼彈阻力值，作出圖7：在同一馬赫數(shù)下，尾翼彈阻力隨著海拔高度的增加在減小，同一海拔高度處，尾翼彈阻力隨著馬赫數(shù)的增加近似線性增加。

模糊控制、專家控制、神經(jīng)網(wǎng)絡(luò)控制、PID控制、自適應(yīng)控制等是工業(yè)機(jī)器人運(yùn)動(dòng)控制系統(tǒng)的主要方法，為了達(dá)到更好的動(dòng)態(tài)性能，最終實(shí)現(xiàn)控制要求，一些新的控制方法就出現(xiàn)了。在控制算法方面既要保證一定的魯棒性，又要滿足參數(shù)實(shí)時(shí)變化的要求。

圖7 尾翼彈阻力隨海拔高度的變化曲線

圖8 阻力系數(shù)隨海拔高度的變化曲線

圖9 阻力系數(shù)隨馬赫數(shù)的變化曲線

在Ma=0.6到Ma=3全聲速范圍內(nèi)，阻力系數(shù)隨著海拔高度的升高而增加，在同一海拔高度處，阻力系數(shù)隨著馬赫數(shù)的變化趨勢(shì)相同。

由數(shù)值方法計(jì)算得到彈體各部分阻力系數(shù)變化值，作出圖11所示曲線，并可以得到：隨著海拔高度的增加，尾翼彈頭部，彈體部，尾翼部的阻力系數(shù)都隨著海拔高度的增加而增加，其中尾翼部的阻力系數(shù)對(duì)海拔高度最為敏感，而彈底阻力系數(shù)隨海拔高度增加而降低。因此尾翼部的阻力系數(shù)升高是造成阻力系數(shù)隨著高度增加的主要原因。

圖10 不同海拔高度處的速度云圖

圖11 Ma=0.6時(shí)彈體各部分的阻力系數(shù)變化曲線

從圖8曲線得到：在海拔高度17 km內(nèi)，阻力系數(shù)與海拔高度基本為線性關(guān)系；亞音速范圍內(nèi)的阻力系數(shù)增加速率要高于超音速范圍的增加速率；在海拔高度超過17 km后，亞音速段，跨音速段的阻力系數(shù)增加較快，而超音速段的增長(zhǎng)速率基本不變，亞音速段，跨音速段的阻力系數(shù)對(duì)海拔高度的變化更為靈敏。

3 阻力系數(shù)擬合

文獻(xiàn)[17]提出：用Logistic曲線與三次拋物線分倆段擬合亞音速段(1.1Ma為界)，超音速段彈丸阻力系數(shù)，該方法擬合誤差較低且分段數(shù)較少，方便于外彈道編程使用。

擬合亞音速段Logistic曲線解析表達(dá)式為：

式中：a分別為L(zhǎng)ogistic曲線下漸進(jìn)線縱坐標(biāo)值；b為上下兩條平行的漸進(jìn)線之間的距離；c，d為擬合曲線待定系數(shù)；x為彈丸飛行馬赫數(shù)；y為彈丸零升阻力系數(shù)。

擬合超音速段三次拋物曲線解析表達(dá)式為：

y=Ax3+Bx2+Cx+D

式中：x為彈丸飛行馬赫數(shù)；y為彈丸零升阻力系數(shù)；A、B、C、D為擬合曲線待定系數(shù)。

用確定系數(shù)R-square來表征擬合的程度：

由上式可知R-square越接近1，表明所擬合方程的擬合程度高。用Matlab擬合方程R-square值均為0.99以上，擬合各海拔高度處阻力系數(shù)方程如表1，擬合曲線如圖12。

圖12 海拔0 km、10 km、20 km處擬合方程曲線

表1 各海拔高度處阻力系數(shù)擬合方程

4 外彈道計(jì)算

射角θ0=10°；初速v0=900m/s；t=0,x=0；y=0km,10km,17km; vx=vvcosθ0；vy=v0sinθ0

由質(zhì)點(diǎn)外彈道計(jì)算得到各海拔高度處發(fā)射該尾翼彈的外彈道特性曲線，質(zhì)點(diǎn)外彈道方程組[18]：

由質(zhì)點(diǎn)外彈道方程計(jì)算得到彈道高，彈丸存速，彈道傾角隨海拔高度的變化值，作出圖13。

圖13 不同海拔高度處彈道諸元的變化曲線

由圖13可以得到：在不同海拔高度處，同一射角，初速下發(fā)射的彈丸彈道軌跡，彈丸速度，彈道傾角有較大的變化。因阻力值隨著海拔高度的增加而降低，所以彈丸存速隨著海拔高度的升高而增加，彈道高隨著海拔高度的增加而下降，彈道傾角隨海拔高度的增加而下降。

5 結(jié)論

在全聲速范圍內(nèi)，同一馬赫處，阻力系數(shù)隨著海拔高度的增加而增大，且基本為線性關(guān)系；尾翼部的阻力系數(shù)隨著海拔高度升高是造成整彈阻力系數(shù)升高的主要原因；在同一海拔高度處，阻力系數(shù)隨著馬赫數(shù)的變化規(guī)律相同。

全聲速段各海拔高度處阻力系數(shù)的擬合方程具有較高的擬合精度。

在相同射角，相同初速，不同海拔處發(fā)射的尾翼彈的彈道軌跡，速度，彈道傾角不同，因其阻力值隨著海拔高度的增加而降低，彈丸存速隨著海拔高度的增加而增加，彈道高，彈道傾角隨著海拔高度的增加而下降，需考慮海拔高度對(duì)彈道參數(shù)的影響。