旋轉彈錐形運動氣動特性研究*

龐川博,康 順,趙 超,蔣勝矩,黨明利

(1 西安現代控制技術研究所, 西安 710065； 2 西安飛行自動控制研究所, 西安 710065)

0 引言

彈箭旋轉飛行具有提高穩定性、消除推力偏心和質量偏心等優點,并且其控制系統更加簡化[1]。在旋轉過程中,彈體周圍流場相對其合攻角平面不對稱,在合攻角平面垂直方向上將產生額外馬格努斯力與力矩[2](也稱面外力與面外力矩)。馬格努斯效應除了引入偏航方向上的氣動干擾外,還會誘發旋進與章動現象,在上述運動耦合作用下,可能導致彈箭發生馬格努斯不穩定而失穩。國內外火箭彈飛行試驗中多次出現因馬格努斯不穩定而失敗的例子,因此旋轉彈箭的氣動特性也成為國內外學者、工程技術人員研究的焦點。準確預測彈體旋轉繞流場的特性，了解其流動機理成為了旋轉彈箭氣動外形設計、控制系統設計的必然需求。

目前旋轉彈箭氣動特性的分析手段可分為風洞試驗、飛行試驗與非定常數值計算。由Dupuis和Hathaway[3]開展的一系列自由飛試驗取得了較好的效果,但試驗中難以獲取流場的細節,無法進一步對旋轉運動中產生的流動現象進行分析,且試驗還具有成本高、周期長等缺點。隨著計算機硬件的發展,非定常數值方法應用于動態過程模擬成為可能,近年來國內外圍繞旋轉彈氣動特性已開展大量工作,其中2007年DeSpirio[4]采用CFD++對M910旋轉彈丸繞流場進行數值模擬,結果表明在亞跨音速條件下RANS/LES混合方法能取得與試驗值符合較好的結果；2012年Vishal[5]采用CFD++中k-ε湍流模型開展了非定常數值計算,計算結果除軸向力偏高外,其余氣動力與動導數都與實驗結果符合較好。在國內,薛幫猛[6]等通過對SOCBT彈丸旋轉開展數值研究,發現尾部形狀對馬格努斯力有較大影響,且馬格努斯力矩與轉速呈線性關系。雷娟棉[7]等通過滑移網格技術計算了高速彈丸的繞流場,發現船尾處邊界層增厚,在船尾處產生的壓力差主要構成了馬格努斯力。肖中云[8]等通過非定常數值計算重點比較了標準形狀和船尾形狀兩種布局的馬格努斯力和力矩的變化可以歸結為亞聲速的繞拐角流動和超聲速的膨脹加速流動對尾部邊界層位移厚度的影響。在2017年石磊[9-10]等對有翼導彈高速旋轉繞流場進行數值模擬,發現大攻角時舵面馬格努斯力占較大比重,同時也針對不同湍流模型對帶翼旋轉彈箭開展了計算,結果表明S-A與SST模型都能較好模擬馬格努斯效應。目前國內外大多研究僅圍繞彈體自轉展開,對彈體旋轉引起的章動、錐動等現象進行數值計算分析的工作還較少,文中在建立描述彈箭多軸角運動方法的基礎上,采用非定常N-S方法結合RANS湍流模型,對標模Basic finner錐動過程開展動態數值模擬,對旋轉過程中流動現象進行分析。

1 數值模擬方法

1.1 控制方程

(1)

式中：Ω為控制體,?Ω為控制體微元面積,n為積分面的單位法向矢量；Q為守恒形式的狀態變量,F(Q)為對流項通量,G(Q)為粘性項通量。在流動未出現較大分離時,非定常RANS方法依然適用于求解旋轉彈箭繞流流場[11],故文中采用Menterk-ωSST湍流模型,其具體構造與特性見文獻[12]。

1.2 多軸旋轉運動的數學描述

根據楊樹興[13]等人總結的旋轉彈箭理想情況下,由初始擾動和陀螺效應引起的二圓運動的復攻角方程,將其成復攻角形式并引入陀螺穩定因子后,得到表達式：

(2)

Δ=C1er1t+C2er2t

(3)

式中r1與r2為特征根,具體寫為：

(4)

利用歐拉公式,式(2)可在α-β平面內寫成隨時間變化的跡線形式,考慮到C1與C2是由初始條件決定的待定常數,將初始條件和運動過程寫成α-β投影平面內幅值與相位關聯的形式,即

Cj=Kjeiφj 0，erjt=e-iφ′jt(j=1,2)

(5)

式(2)最終寫作

Δ=K1eiφ1+K2eiφ2

(6)

由式(6)可以看出,在復平面內,通解由兩個線性無關的不同幅值、不同頻率的圓運動組合而成,依照頻率的高低,分別稱之為快圓運動與慢圓運動。考慮彈體自轉時,二圓運動實際上為彈箭繞3個不同轉軸旋轉的復合角運動。選取合適的參考系對不同角運動進行描述,再通過坐標變換,即可獲得多軸角運動的描述形式。

在相應參考系內,定義沿周向旋進過程掃過的角度為進動角,沿徑向臂長變化的角度為章動角,上述兩個角度及其微分形式即可對單個圓運動進行完整描述。描述二圓運動所需角度參數包括自轉角φ0(對應自轉運動)、快圓進動角φf與快圓章動角Kf(對應快圓運動)、慢圓進動角φs與慢圓章動角Ks(對應慢圓運動),下面建立相應參考系對各個角運動進行定義。

圖1 慣性坐標系以及初始彈體姿態時尾翼編號

彈體系與慣性系是彈道中常用參考系,其中彈體系與彈體固連,慣性系與地面固連。令慣性系xo方向與來流方向保持一致,且初始時刻彈體軸線處于慣性系xoy平面內(見圖 1(a)),此時對4片尾翼進行編號(見圖 1(b))。建立與彈軸固連但不跟隨其旋轉的準彈體系,得到用于描述彈體自轉的自轉角φ0；將準彈體系繞z軸旋轉Kf后得到快圓坐標系,在此基礎上類似建立準快圓坐標系,得到用于描述快圓旋進運動的進動角φf；將準快圓坐標系繞z軸旋轉Ks后得到慢圓坐標系,同理得到在準慢圓坐標系中用于描述慢圓旋進運動的進動角φs。易知最終的準慢圓坐標系與慣性系是重合的。

注意上述各角度均在各自其給定的參考系下進行描述,目的是使其物理意義更清晰。將自轉運動、慢圓運動、快圓運動在各自最初定義的坐標系中分別進行坐標轉換(其中自轉運動需將坐標系旋轉4次,快圓運動需將坐標系旋轉2次,慢圓運動旋轉0次),最終可寫出慣性系下復合角運動各分量的表達式：

(7)

其中：

式中：φf0與φs0為快圓運動與慢圓運動在初始時刻所處姿態角,當其為零時彈體軸線位于慣性系xoy平面內且彈尾位于y+方向；w0、ws與wf分別為自轉角速度、慢圓運動和快圓運動的旋進角速度。式(7)具有較好的普遍性,適用于描述定點三軸及以下多軸角運動的姿態及運動過程。

在進行數值計算時,將角運動分量投影至不同時刻彈體坐標系中往往更有實際意義。首先引入角運動的瞬軸概念,其物理意義為：剛體繞某一矢量進行角運動時,若該矢量方向決定了剛體轉動的方向,該矢量大小決定了剛體轉動的角度,則定義該矢量為當前時刻剛體轉動的瞬軸。瞬軸與角位移滿足關系：單位時間角位移=瞬軸×單位時間。

容易知道慣性系下瞬軸可寫為：

(8)

式(8)為了形式簡潔,只將角速度寫成時間相關變量。上標I表示該變量在慣性系下進行描述,B1為t1時刻的彈體系。根據上文中已得到的式(7),目前可以寫出在慣性系下每個時刻的瞬軸。當彈體旋轉時,物面及空間網格一同跟隨旋轉,旋轉后得到新的參考系,稱為當地隨體坐標系。

在初始時刻,慣性系與當地隨體坐標系重合,此時瞬軸在兩個坐標系中表達式相同,即

sB1(t)=sI(t)

(9)

下一時刻,網格將繞瞬軸進行第一次旋轉,由于瞬軸已知、單位時間步長已知,可構造第一個旋轉矩陣,記為R(t1)。旋轉矩陣通過矩陣運算將目前坐標系(舊時刻的當地隨體系,在初始時刻與慣性系重合)旋轉至新坐標系(下一時刻的當地隨體系)。

容易得到下一時刻瞬軸在本時刻當地隨體系中表達式

sB2(t2)=R-1(t1)·sB1(t2)=R-1(t1)·sI(t2)

(10)

以此類推,再下一時刻

(11)

最終可以寫出慣性系下瞬軸旋轉至任意時刻當前體軸系下的表達式:

(12)

此時便將剛體的旋轉運動方程寫成了離散的、投影在當前時刻當地隨體參考系中的表達形式,適用于非定常動態數值計算中運動過程的描述。

1.3 計算模型

計算模型采用Basic Finner導彈標模,基本外形為十字布局的帶翼導彈,由彈頭錐段、圓柱段與尾翼組成,質心位于距彈頭55%處,基本幾何尺寸見圖 2。該標模為美國空軍在二戰后為開展無控帶翼導彈空氣動力學和飛行力學特性研究而設計的模型,目前已經過大量的靜態與動態氣動試驗研究,有豐富的靜態和動態試驗及計算數據。

圖2 Basic Finner基本外形尺寸

1.4 數值方法驗證

共采用3套疏密不同的網格用于考察數值方法的計算準確性,分別為1)粗網格：網格總量大約400萬,邊界層底層高度為5×10-6m；2)中等網格：網格總量大約600萬,邊界層底層高度為2×10-6m；3)較密網格：網格總量大約800萬,邊界層底層高度為2×10-6m。計算網格采用多塊對接的混合網格,彈身表面、近壁面周圍以及尾跡區域采用加密后的六面體網格以保證流動特征的捕捉,遠場區域采用四面體網格以節省計算資源。計算網格如圖 3所示。

圖3 計算網格示意圖

計算馬赫數取1.254,不同網格計算結果均與試驗值[14]進行對比。結果如表1所示,表中第一行從左到右分別為：零攻角條件下法向力系數斜率、軸向力系數、相對質心的俯仰力矩系數對攻角的斜率以及滾轉力矩系數。所有計算均收斂較好,加密網格后,中等規模網格計算值與試驗值已非常接近,可認為中等規模網格已滿足需求。

表1 不同疏密網格計算結果對比

2 計算結果與討論

2.1 計算工況說明

由式(7)可知,當Kf與φ′f為零時,二圓運動退化表現為錐形運動,此時慢圓運動轉軸為錐動軸,慢圓運動的章動角為錐動角,慢圓運動的進動角速度即為錐形運動的旋進角速度；文中主要計算狀態為錐動角Ks=10°,不同自轉角速度與進動角速度下彈體的錐形運動,具體見表2。

表2 不同工況彈體運動參數

工況1、2、3的進動角速度方向和尾翼斜置產生的滾轉力矩方向一致；工況1、2彈體無自轉,此時運動方式為似月運動,彈體表面迎風區域與背風區域保持不變,相對來流的姿態也不變。工況3中彈體繞自身彈軸旋轉,自轉角速度是進動角速度的兩倍,此時彈體相對來流的姿態在“十”和“X”之間不斷變化。非定常時間步長Δt=3.87×10-5s,進動一周的周期T=0.013 9 s,非定常時間推進360步。自由來流狀態為：馬赫數Ma=1.254,雷諾數Re=8.76×106,溫度Tref=288.15 K。

為方便說明,后文將基于慣性系與隨體風軸系進行分析,兩個坐標系定義如圖4。

圖4中左圖為來流與彈體初始姿態關系,右圖為兩個坐標系的相對關系。原點o位于彈體質心處,z軸滿足右手定則,xoy不隨彈體運動而改變,來流與ox方向平行；彈軸與ox夾角為合攻角,由慢圓運動章動角和快圓運動章動角組成(上文中已規定其在初始時刻位置)。彈體作錐形運動旋進時,設坐標系xbybzb跟隨彈體進動而旋轉(并不跟隨自轉而旋轉),且obxb與ox始終重合。

圖4 慣性系xoy與隨體風軸系xbybzb示意

2.2 彈箭似月運動氣動特性討論

在上文給出的慣性系下將彈體受力分別沿oy方向和oz方向進行分解,得到的Cy與Cz隨進動角變化的曲線如圖 5所示。可以看出彈體在不同位置受力滿足正弦變化規律,不同的周期之間結果也保持良好一致,且兩個工況計算結果相差較小。

圖5 慣性系下Cy與Cz隨不同相位角變化

分別取上述兩個工況非定常計算的最后一個周期,與不同相位下的定常結果以及對非定常結果進行擬合的曲線進行比較,兩個工況有類似的結果,故此處僅以工況2為例,給出曲線變化如圖 6所示。擬合曲線與非定常計算結果符合較好,可用于描述非定常曲線的相位；與定常結果相比,動態過程受力最大值更大,且具有相位提前的特點,即在進動過程中,在進動角尚未達到179°時(尾翼斜置的影響),Cy已達到正向最大值,Cz同理。通過對比擬合曲線和定常曲線可以得到相位提前的角度,工況1的法向分量Cy的相位提前角大約為5.296°,側向分量Cz的相位提前角大約為5.094°；工況2的法向分量Cy的相位提前角大約為6.724°,側向分量Cz的相位提前角大約為6.68°。

相位提前現象與給定的進動角速度有關,下面在隨體風軸系下進行討論,其中Cyb為合攻角平面內的法向分量,Czb為垂直于合攻角平面的側向分量。彈體繞過質心的錐動軸(即慢圓進動軸)轉動構成了一個旋轉系統,在不對稱斜置尾翼的作用下,該旋轉系統在當前馬赫數下具有固定的平衡轉速,當進動角速度等于當前旋轉系統的平衡轉速時,滾轉力矩Mxb為零,即Czb與Cyb在滾轉方向產生的力矩相互抵消,已知Cyb主要由來流在合攻角方向產生,其數值遠大于Czb,可認為此時Czb對合力方向影響很小,幾乎不發生相位提前效應；當進動角速度大于當前旋轉系統平衡轉速時,由于滾轉力矩Mxb不等于零,彈體Czb必然不為零,此時Czb與Cyb在當前位置合力的方向與之后某時刻所處相位的法向力方向一致,故出現相位提前的特征(見圖 7所示)。進動角速度越大,相應的面外力也越大,相位提前效應就越明顯,非定常計算結果合力的最大值也越大(合力由合攻角平面內的Cyb與面外力Czb所合成)。

圖6 Cy與Cx相對定常解的相位提前量

圖7 似月運動引起的相位提前現象

從上文中可知,Cy與Cz曲線的相位提前角并不完全相同,兩條曲線仍具有一定相位差,Cyb與Czb將呈小幅振蕩趨勢(若Cy與Cz始終同相位,則Cyb與Czb為固定常值)。將兩工況非定常解投影至隨體風軸系下,得到如圖 8所示結果,分析圖中結果得到以下結論：在當前所采用的網格分布和湍流模型的條件下,彈體沿隨體風軸系分解得的法向力Cyb和側向力Czb均隨進動角變化產生小幅振蕩,進動每經過90°相位,Cyb與Czb完成一個振蕩周期；進動角速度變化對Cyb與Czb周期和振幅影響較小,但對Czb的均值有一定影響,當進動角速度大于旋轉系統平衡轉速時,Czb均值隨進動角速度增大而增大。

在圖9觀察工況2在相位角為零時尾翼1、3流場分布,由于尾翼斜置的原因,尾翼1的進動背風面同時也是自由來流相對尾翼斜置方向的迎風面,尾翼3的進動背風面則是來流相對斜置方向的背風面,且尾翼1翼梢處線速度更大,在來流與旋轉綜合作用下,尾翼1背風區域能觀測到流動分離現象,產生的渦附著于翼梢進動方向背風面,來流的影響使得越接近后緣,渦規模發展越大,分離趨勢更加明顯,而尾翼3翼梢處并沒有明顯的分離現象。尾翼1、3的受力變化很大程度上決定了Cyb與Czb的變化趨勢,但由于振蕩幅值和渦的尺寸均很小,進一步觀測和分析渦發展和脫落需要更精細的網格及采用LES方法。

圖8 隨體風軸系下受力結果

圖9 零相位時尾翼附近流場渦量圖

2.3 錐形運動中自轉的影響

工況3考慮自轉角速度,彈體旋進過程中自身相對姿態也在不斷變化,取彈體底部過彈軸一點與翼梢處某固定點作為觀察點,彈體進動一個周期的觀察點軌跡如圖 10所示(為方便觀察,示意圖中旋轉中心設在彈頭,實際計算時彈體旋轉中心在質心處),紅色虛線為進動過程的軌跡,而藍色虛線則為彈體繞彈軸自轉的軌跡,在一個進動周期內彈體自轉兩圈。圖 11為工況3與工況2全彈在慣性系下受力分量隨進動角變化曲線,通過對比,工況3曲線仍具有一定的正弦變化特征,但與似月運動結果已有較大差別：1)相位略有提前；2)在曲線的波峰與波谷處,即相位角90°、180°與270°處,工況3與工況2的姿態完全相同,但受力分量明顯較小；3)曲線整體光順性較差。

圖10 彈體旋進一周期尾部與翼梢觀察點軌跡

圖11 工況2、3在慣性系下受力曲線

在隨體風軸系下進一步分析,首先給出工況3與工況2的受力曲線如圖 12所示。對比可知,彈體繞彈軸自轉后,其合攻角平面內Cyb與合攻角平面外Czb的均值、振幅、周期均有較大改變。對于工況3,每進動180°,曲線經過一個大周期,每個大周期內有4個小周期,每個小周期對應彈體自轉90°；工況3中受力曲線振幅遠大于工況2的似月運動,說明彈體自轉誘發的翼梢處渦脫落現象更加顯著,翼梢處的流動分離也將引起合攻角方向上升力的損失,因旋轉引起的平均升力損失大約占總升力的4.55%；兩種工況側向力Czb的變化主要由彈體自轉產生的馬格努斯效應造成,在有翼彈箭翼體干擾和彈體進動過程與氣流的相對作用綜合影響下,側向力Czb均值明顯增大,增量占全彈法向力的6.5%左右。Cyb與Czb合力的變化與振幅的不規則改變可以看作是馬格努斯效應影響的結果。

圖12 隨體風軸系下工況2與工況3受力曲線

圖13給出進動90°(即1/4周期)時不同截面的渦量分布。兩種工況彈身背風區均出現明顯流動分離現象；在彈身圓柱段,似月運動使背風區下側分離渦受到擠壓更貼近壁面,而上側則受相對似月運動的等效來流影響,分離渦遠離彈體壁面。工況3結果見圖 13(b),彈體旋轉引起的環流在彈體上側與來流分量相反,二者相互阻礙使分離渦更遠離彈體且更容易脫落；彈體下側環流與來流方向相同,使邊界層局部速度梯度更大。

圖13 t=1/4T時刻不同截面渦量分布

在x=0.54L截面作工況3的壓力分布云圖,同時給出工況2、3在該相位下彈身截面處的流向壓力分布曲線,如圖 14所示。從圖中可知,由于彈體的進動,兩個工況彈身下側壓力均大于上側壓力,同時在背風面處都表現出了明顯的上下不對稱現象,最高壓力點從270°處偏移至255°處；工況3中由于彈身的自轉影響,背風區壓力分布不對稱性更加明顯,且普遍在180°相位至270°相位壓力更大,彈身旋轉所產生的實際面外力指向圖中40°相位方向處,符合理論上旋轉彈彈身產生馬格努斯力方向的規律。靠近彈尾處流動分離與不對稱性更加顯著,由于彈翼的遮擋作用,背風面漩渦強度降低,渦的形狀分布特點與彈身類似。翼梢渦的規模取決于當地垂直翼面氣流的相對速度,在尾翼中部靠中后緣區域受激波干擾,邊界層增厚,流動更易分離(見圖 15)。

圖14 x=0.54L壓力分布云圖與曲線圖

圖15 尾翼不同截面馬赫數分布

3 結論

文中針對旋轉誘發的復雜角運動現象,建立適用于計算流體力學的多軸角運動表達形式,對標模Basic Finner的似月運動以及錐形運動開展數值計算,從流場結構和氣動載荷變化進行了研究,通過分析得出以下結論：

1)彈箭作似月運動時,沿慣性系各方向受力隨相位角變化的曲線,相比定常結果具有相位提前的特征,且進動角速度越大,相位提前量越大。

2)錐動過程中尾翼翼梢及彈尾處存在小尺度流動分離,使得彈體在合攻角平面內與平面外受力呈振蕩現象；彈體自轉速度為零時側向力振蕩幅值較小,當自轉速度不為零時,側向力均值與振蕩幅值顯著增大。進一步準確觀測和分析該現象需要更精細的網格與更適合的湍流模型。

3)彈體自轉會使得彈身背風面以及尾翼翼梢處不對稱流動分離現象更顯著,二者綜合作用下彈體產生額外的側向力。