基于計算流體力學與剛體動力學耦合的高速旋轉彈丸彈道計算方法

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

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

0 引言

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

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

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

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

1 CFD數值計算方法

1.1 基于彈軸運動的ALE形式控制方程

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

(1)

(2)

Fc為靜態網格下的對流通量，v?Ω為附加在控制體表面的反變速度。

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

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

(3)

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

1.2 修正簡單低耗散迎風矢通量分裂格式

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

(4)

(5)

cij=(cl+cr)/2,

(6)

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

壓力通量為

(7)

(8)

(9)

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

1.3 邊界條件、時間推進及湍流模型

(10)

(11)

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

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

2 基于CFD的彈丸飛行動力學模型

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

Lac=diag[-1 1 -1].

(12)

第1彈軸坐標系到彈道坐標系的轉換矩陣為

(13)

(14)

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

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

(15)

彈丸質心運動的運動學方程組[1]為

(16)

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

(17)

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

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

(18)

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

(19)

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

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

(20)

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

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

(21)

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

彈道方程組中繞質心轉動的運動學方程組和δa、δd、β3個角關系方程參見文獻[1]，不再贅述。

3 CFD/RBD耦合計算方法

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

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

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

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

(22)

式中：

(23)

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

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

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

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

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

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

圖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中采用雙時間步推進到的tn+h時刻為n+1時間層面。

4 仿真驗證與分析

4.1 M549彈丸定態運動和規定運動模擬

計算模型采用美國M549旋成體彈丸，具體結構尺寸參見文獻[24]。利用Solidworks商業軟件進行幾何建模，如圖3所示。

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

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

圖4 M549彈丸計算域網格Fig.4 Computational domain mesh of M549 projectile

通過使彈丸做定態飛行和規定運動獲得彈丸的氣動特性，將結果與Spinner數據[24]進行對比，以檢驗本文CFD數值方法的可靠性。

(24)

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

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

總之，在馬赫數為2.0時，靜態氣動系數和動態氣動系數的仿真結果相對Spinner數據的誤差分別小于6%和10%，數值仿真的氣動系數較可靠。

4.2 M549彈丸自由飛行仿真驗證與分析

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

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

對于傳統6自由度彈道計算，時間步長通常取5 ms左右[1]。為了研究3種力和力矩耦合方法求解彈道的能力，取步長h=5 ms進行仿真，得到攻角關系曲線如圖6所示。通過觀察圖6中的攻角變化情況發現：采用定值法效果最差，仿真攻角結果發散，計算彈道失?。徊逯捣ê途o耦合法均模擬出一般旋轉穩定彈丸的雙圓運動，其中緊耦合法計算結果與基準結果更加接近，即緊耦合精度最好。

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

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

圖7 彈道參數仿真結果Fig.7 Simulated results of trajectory parameters

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

4.3 修正組件不轉時二維彈道修正彈模擬與分析

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

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

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

圖9 攻角仿真結果Fig.9 Simulated result of angle of attack

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

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

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

5 結論

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

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

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

3)緊耦合方法模擬具有復雜外形的旋轉穩定彈飛行時，能夠反映出彈丸的運動特性，流場變化規律符合客觀實際。

本文提出的耦合計算方法可為新型旋轉穩定彈，如加裝了精確制導組件的二維彈道修正彈設計和計算提供一種新的理論和方法。今后可進一步考慮控制部件的運動，為新型旋轉穩定彈控制組件的設計提供基礎。