低速旋轉導彈滾轉特性數值模擬研究

徐科杰，王學德，張 超，張雋研，王鋅晨

(南京理工大學 能源與動力工程學院， 南京 210094)

現代制導武器中滾轉控制是非常重要的問題，無論是帶有控制系統的戰術導彈、火箭彈、炮彈都需要控制其滾轉速度，以此測算出飛行滾轉角達到導彈飛行姿態控制的目的，從而提高導彈飛行穩定性和打擊精度。導彈的穩定轉速(即平衡轉速)與導轉力矩以及滾轉阻尼力矩息息相關，因此研究旋轉導彈的滾轉特性極其重要[1]。

從20世紀70年代開始，國外學者利用風洞試驗對導彈滾轉特性進行了大量研究。F.J.Regan、Leroy、Collins.J.A[2-4]等對標準模型進行了風洞實驗并得到了一些有用結論。隨著計算機性能的大幅度提高，研究人員開始采用數值模擬方法對滾轉特性進行研究。Pual Weinacht、Ameer G.Mikhail等[2-4]驗證了基于N-S方程的數值計算方法計算導彈氣動特性。Erdal Oktay、John F.Stalnaker等[6,7]采用基于不同算法的求解器對標準模型Basic Finner進行了非定常計算。在國內，鄧帆[8]發現在超聲速條件下的氣流雍塞是導致柵格翼導彈的滾轉阻尼導數減小的主要原因。周德娟[1]分別采用風洞試驗和數值仿真對不同翼身組合導彈滾轉特性進行了系統研究，兩者計算結果吻合較好。余奇華、敬代勇[9]對不同轉速下的旋轉尾翼鴨式布局彈箭的流場特性進行了數值模擬，發現尾翼旋轉可以顯著提高其滾轉控制能力。

國內外學者從風洞試驗與數值仿真兩方面對導彈的滾轉特性進行了大量研究，基本上以研究馬赫速、攻角與導彈滾轉特性之間的關系為主，但在滾轉特性隨轉速以及尾翼斜置角的變化方面研究較少。本文基于滑移網格技術，利用計算流體力學(CFD)軟件對不同尾翼斜置角以及轉速下導彈三維繞流流場進行了數值模擬，從流場結構以及氣動特性兩部分對導彈的滾轉特性進行了闡述。通過數值模擬計算不同轉速下的滾轉力矩系數，使用Matlab進行數據擬合得到了滾轉力矩系數隨轉速變化的一次擬合函數，進而得到導彈在不同工況下的平衡轉速，并且與工程算法得到的平衡轉速進行對比，得到快速計算導彈的滾轉力矩系數以及平衡轉速的經驗公式。

1 數值模擬方法

滑移網格技術是20世紀80年代逐漸興起的一種基于非穩態思想的數值模擬方法，近年來得到了廣泛的應用，尤其是在彈箭旋轉的氣動仿真領域。圖1為滑移網格技術的信息交換示意圖。整個流動區域被分為包含1、2、3單元和4、5、6單元的兩個不同的區域，其中一個區域相對另一個區域做相對滑移運動來模擬真實的流動過程。ABCD和EFGH是各自子區域與滑移面重合的邊界面。在數值模擬過程中，1、2、3單元和4、5、6單元在運動過程中，會形成abcdefg節點，并沿著交界面做相對滑移運動。其中子塊2的信息通過面cd和面de傳遞到子塊4和子塊5，其他單元之間的信息傳遞與此類似。

滑移網格技術作為一種特殊的動網格，適用于周期性旋轉運動物體的流場計算，具有計算精度高，占用內存少等優勢[10]。當模擬彈箭旋轉運動時，只需要一個非旋轉區域和一個包圍彈體的旋轉區域，兩者之間存在一個交界面，交界面的網格節點數不必相同。僅僅在滑移網格交界面進行數據插值，即可保證兩個區域之間通量守恒。圖2(a)和圖2(b)給出了網格區域劃分示意圖。

采用基于滑移網格下的三維積分守恒型非定常N-S方程組[11]，由質量守恒方程、動量守恒方程和能量守恒方程組成。

式中：ρ,U,e分別表示氣流的密度、速度和單位體積的能量；vm為彈體的運動速度；P為靜壓力；I為單位張量；τ為黏性張量；q為熱通量；S和v分別表示控制體的體積分和面積分的積分區域；r為S的外法向單位向量。

采用有限體積法對空間進行離散，對流項采用二階通量差分分裂(FDS) AUSM+格式，黏性項采用中心差分格式。

Realizablek-ε模型是對RNGk-ε湍流模型做了修正而產生的，因此可以更為準確地反映流體真實的流動過程，尤其是旋轉流動計算中氣流的流動狀態。

2 算法有效性驗證

2.1 驗證模型

選取美國陸-海軍動導數計算標準模型(ANF)進行數值算法的有效性驗證，ANF模型經過大量風洞實驗和飛行試驗驗證，被廣泛采用為動導數的計算標準模型，如圖3所示，ANF-1與ANF-2模型的差別是：ANF-1無尾翼斜置角而ANF-2有尾翼斜置角。表1和表2給出計算條件，其中攻角(α)、馬赫數(Ma)、雷諾數(Re)、總壓(P0)、靜壓(P)和總溫(T0)、靜溫(T)同風洞實驗條件,Ω·為對應的導彈無量綱旋轉速度(Ω·的定義為Ω·=D*Ω/V,D為彈體直徑,V為當地聲速，Ω為導彈轉速(rad/s),N為外循環中每轉動一周的總計算的迭代次數,Δt為外循環的計算步長,i為每一次外循環計算中內循環的迭代次數。

MaRe(106)P0/PaT0(k)Ω·2.491.8644 816311.10.034 6NΔt(10-6)iD/mm1 4401.1312031.75

表2 ANF-2模型的計算條件

利用ICEM商業軟件，采用多塊對接網格生成技術對整個計算區域生成六面體網格。縱向采用C型網格，橫向采用O型網格。第一層距離物面為1.0×10-5m，網格伸展比為1.1，控制第一層網格到物面的距離和網格伸展比的目的是為了保證y+<1。圖4(a)和圖4(b)分別是ANF-1模型彈體頭部和尾部附近網格。

彈箭表面設置為無滑移壁面條件，彈體跟隨旋轉區域以相同的轉速同步運動。非旋轉區域和旋轉區域的交界面設置為滑移邊界條件；非旋轉區域的外邊界設置為壓力遠場條件。

2.2 算法驗證分析

圖5(a)和圖5(b)分別給出了不同攻角時ANF-1模型的滾轉力矩動導數和俯仰力矩系數的計算值與實驗值[12]。從圖5可以看出滾轉力矩動導數和俯仰力矩系數與相應的實驗值吻合較好，最大誤差分別為17.6%和8.62%(最大誤差均出現在攻角α=80°附近)。同時發現攻角越大，誤差越大。這是由于攻角越大，自旋導彈附近的流場變得更為復雜。為了進一步驗證該數值模擬方法的精確性，對比了ANF-2模型的軸向力系數和法向力系數隨馬赫數的變化規律，并且與實驗值進行了比較，如圖6所示軸向力系數與法向力系數與實驗值[13]吻合較好，其最大誤差分別為3.3%(1.1馬赫數下)和7.2%(3.5馬赫數下)。可知該數值模擬方法在求解小攻角、全馬赫數情況下有較高的精確度，驗證了算法的正確性。

3 某低速尾翼彈滾轉特性計算

計算模型為某尾翼彈，由錐形部、圓柱段和尾翼組成，其尺寸外形如圖7所示,彈徑為0.4 m，彈長為L為7.3D。選取表3中的參數作為計算狀態，氣動力和力矩計算的參考點位于彈體質心位置，參考面積為彈體最大橫截面積，參考長度為彈徑D。

3.1 網格尺度以及時間步長獨立性驗證

模型網格數目以及物理時間步長的選擇對仿真計算結果有較大影響，尤其是在非定常計算中。若網格數目較少，將無法準確地捕獲到流場信息甚至導致計算出錯；若網格數目過大，則消耗計算資源，延長計算周期。與網格獨立性原理一樣，如物理時間步長太大，則會導致計算發散，如物理時間步長過小的話，又會加大計算時間。因此需要選取一個合理的時間步長以及網格數目。表4給出了Ma=2.0、α=2°、Ω·=0.015 1、γ=0°下，不同網格數下阻力系數CD、升力系數CL、俯仰力矩系數Cm以及滾轉力矩系數Cl參數的計算結果。可見隨著網格的不斷加密，該導彈的氣動特性參數逐漸趨于收斂。考慮到900萬與800萬網格下的計算結果差異較小(以900萬網格下的氣動數據作為基準，四組氣動力參數的最小差值分別為1.3%、1%、1.2%、3.8%)，因此后續計算以800萬網格為準。圖8給出了導彈四組氣動特性參數隨物理時間步長Δt的變化曲線(計算條件參考表3所列的計算參數) ，考慮到Δt=0.000 4 和Δt=0.000 6中各氣動特性參數差異較小(以Δt=0.000 4下的計算結果作為基準，四組氣動力參數的最小差值分別為0.5%、1.1%、1.6%、1.3%)以及需要消耗的計算周期，本文采用Δt=0.000 6。

表3 模型的計算參數

表4 網格數對氣動特性參數的影響

3.2 流場分析

如圖9(a)和圖9(b)所示為尾部截面的速度矢量圖，從圖9(a)可知，由于攻角的存在，使得在流場區域存在一個整體向上的速度分量。從圖9(b)可以明顯看出，導彈在旋轉過程中，尾翼引起其附近氣流隨導彈一起順時針旋轉。A區域與C區域氣流速度矢量方向相反，因此形成一個速度矢量為0的區域即B區域，同時也會在翼根位置即D區域形成一個較小的渦流。

3.3 氣動特性參數分析

圖10為馬赫數Ma=2.0、尾翼斜置角γ=1°和不同轉速下的尾翼位置的壓力云圖。從圖10(a)可知，對于上下尾翼，由于整個尾翼處于背風(或者迎風)氣流狀態，當存在較小的尾翼斜置角時，使得上尾翼的兩側翼面再次處于不同的氣流中，產生較大的壓力差。同時發現，上下尾翼產生的壓力差方向相反，形成了加速導彈轉速的滾轉力矩。對于左右尾翼而言，攻角的存在使得其兩側翼面產生較大的壓力差，但由于左右尾翼壓力差方向相同，因此對滾轉力矩貢獻較小。對比圖10(a)和圖10(b)，隨著轉速的增加，上下尾翼的兩側翼面壓力差逐漸變小，由壓力差產生的滾轉力矩變小。從下文的定量分析中也可以得到該結果。

通過數值模擬獲得了自旋尾翼彈在各種飛行條件下的滾轉特性數據，如圖11所示。由圖11(a)可知，滾轉力矩系數的絕對值隨著馬赫數的增大呈現先增大后減小的趨勢，并且在Ma=1.1時其滾轉力矩系數的絕對值最大。并且滾轉力矩系數的絕對值隨著轉速的增大顯著增大；同時發現圖11(a)中的滾轉力矩系數均小于零，這是由于尾翼斜置角γ=0°，導彈無法獲得使其主動滾轉的力矩(即導轉力矩為0)而只存在阻礙其旋轉的滾轉阻尼力矩，故而滾轉力矩系數均為負值。

對比圖11(a)和圖11(b)發現，兩者隨馬赫數的變化趨勢一致：其絕對值均隨馬赫數的增大呈現先增大后減小趨勢。兩者的區別在于：圖11(b)中的滾轉力矩系數均為正值，這是由于尾翼斜置角γ=1°時，該導彈存在一個使其主動滾轉的力矩(即導轉力矩)并且導轉力矩的值大于滾轉阻尼力矩的值，所以其滾轉力矩系數表現為正值，在Ma=1.4附近滾轉力矩系數取得最大值；同時還發現在圖11(b)中，導彈轉速越大，其滾轉力矩系數越小，這是由于轉速越大，導彈的滾轉阻尼力矩也就越大。

對比圖11(b)、圖11(c) 和圖11(d)可知：在同一馬赫速和轉速條件下，尾翼斜置角越大，則滾轉力矩系數越大。這是由于尾翼斜置角越大，尾翼兩側的壓力差也越大，導彈產生的導轉力矩越大。

通過數值模擬獲得了自旋尾翼彈在各種飛行條件下的滾轉特性參數，如圖12給出了不同Ma下導彈滾轉力矩系數Cl隨無量綱轉速Ω·的變化趨勢。

從圖12(a)可知，當尾翼斜置角γ=0°時，導彈滾轉力矩系數的絕對值隨轉速的增大呈現增大趨勢，并與無量綱轉速存在線性關系。對于同一轉速條件下，在Ma=1.1下導彈的滾轉力矩系數(即滾轉阻尼力矩)的絕對值最大，Ma=2.0下導彈的滾轉力矩系數的絕對值最小。同時Ma=1.1和Ma=0.9條件下的滾轉力矩系數曲線基本平行。對比圖12(b)和圖12(c)、圖12(d)可知：其滾轉力矩系數均隨無量綱轉速的增大呈現減小趨勢。

通過Matlab軟件的擬合函數[b,bint,r,rint,stats]=regress(x,X)對導彈滾轉數據進行擬合，得到了不同馬赫數下，滾轉力矩系數Cl與無量綱轉速Ω·以及尾翼斜置角γ的關系式。如表5所示，發現滾轉力矩系數與無量綱轉速Ω·以及尾翼斜置角γ呈線性關系。將擬合值與仿真值對比發現，兩者的最大平均誤差小于9%，因此認為該關系式較為準確地描述滾轉力矩系數與轉速以及斜置角三者之間的關系。

根據表5的擬合函數得到了導彈飛行的平衡轉速(即滾轉力矩系數為零)，同時與工程算法得到的平衡轉速進行對比，結果如表6～表8所示。當尾翼斜置角γ=0°時，旋轉導彈在飛行過程中只存在阻礙其旋轉的滾轉阻尼力矩，即該導彈處于減旋過程中直到轉速為0，因此其平衡轉速為0。如表6所示，當尾翼斜置角γ=1°時，不同馬赫數下的導彈的平衡轉速不同，并且馬赫數越大，其平衡轉速越大。對比表6和表7、表8可知：對于同一馬赫數而言，尾翼斜置角越大，其平衡轉速越大。將數值計算方法得到的平衡轉速和工程算法[14-16]得到的平衡轉速進行對比發現，兩者的最大誤差不超過19%，考慮到所有的工程算法本身具有一定的誤差，因此認為本文中通過擬合得到的經驗公式(即表5中的擬合表達式)真實有效并且可以對平衡轉速進行估算。

表5 擬合表達式

表6 γ=1°時導彈平衡轉速

表7 γ=2°時導彈平衡轉速

表8 γ=4°時導彈平衡轉速

4 結論

1) 隨馬赫數增加，平均滾轉力矩系數絕對值呈現先增加后減小的變化趨勢，與常規滾轉力矩變化規律一致。

2) 當尾翼斜置角為0°時，尾翼不產生導轉力矩，其滾轉力矩系數均為負值即滾轉阻尼力矩；隨轉速增加，滾轉阻尼力矩增強，且滾轉阻尼力矩呈現典型的線性變化規律；當尾翼斜置角為1°、2°、4°時，其滾轉力矩均為正值，且隨轉速的增大而減小。

3) 用Matlab軟件擬合出滾轉力矩系數與無量綱轉速以及尾翼斜置角之間的函數關系式，并且將擬合值與仿真值以及工程算法得到的計算值進行比較，三者吻合較好，表明本文得到的經驗公式真實有效。由于該公式可以快速地求出相應計算條件下的滾轉力矩系數，以及相應的平衡轉速，因此具有較大的實用價值。