尾翼彈高低空氣動特性仿真及外彈道計算

岳 通，王惠源，張成卿

(中北大學 機電工程學院，太原 030051)

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

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

1 數值計算方法

本文研究海拔高度在20 km以內，小口徑尾翼彈的氣動特性，其所屬流域為連續流域，氣體的間斷效應并不顯著，連續性假設成立，所以仍可運Navier-Stokes方程進行計算，選用雷諾平均的N-S方程作為控制方程。

控制方程為[11]：

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

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

尾翼彈外形參數如圖1所示。

圖1 尾翼彈外形參數

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

圖2 計算域網格

圖3 彈體表面網格

由美國1976年標準大氣分段逼近公式[15]計算：

計算位勢高度H,它與幾何高度Z的關系為：

H=Z/(1+Z/R0)

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

1) 當0≤Z≤11.019 km時：

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

2) 當11.0191 km≤Z≤20.063 1 km時：

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

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

3) 黏性系數隨高度變化的計算公式：

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

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

由以上公式計算得到各海拔高度的氣象參數值，作出相關曲線如圖4所示。

圖4 氣象參數變化曲線

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

算法有效性驗證：

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

圖5 標準模型

圖6 標模彈零升阻力系數曲線

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

2 計算結果與分析

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

模糊控制、專家控制、神經網絡控制、PID控制、自適應控制等是工業機器人運動控制系統的主要方法，為了達到更好的動態性能，最終實現控制要求，一些新的控制方法就出現了。在控制算法方面既要保證一定的魯棒性，又要滿足參數實時變化的要求。

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

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

圖9 阻力系數隨馬赫數的變化曲線

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

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

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

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

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

3 阻力系數擬合

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

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

式中：a分別為Logistic曲線下漸進線縱坐標值；b為上下兩條平行的漸進線之間的距離；c，d為擬合曲線待定系數；x為彈丸飛行馬赫數；y為彈丸零升阻力系數。

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

y=Ax3+Bx2+Cx+D

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

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

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

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

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

4 外彈道計算

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

由質點外彈道計算得到各海拔高度處發射該尾翼彈的外彈道特性曲線，質點外彈道方程組[18]：

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

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

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

5 結論

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

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

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