旋轉(zhuǎn)彈體馬格努斯效應(yīng)數(shù)值模擬方法研究

陳白冰，駱振華，袁振宇，蔣曉磊，陳兵

(1.中國(guó)兵器工業(yè)集團(tuán) 西安現(xiàn)代控制技術(shù)研究所，西安 710065) (2.西北工業(yè)大學(xué) 航空學(xué)院，西安 710072)

0 引 言

導(dǎo)彈、火箭彈等武器在飛行過(guò)程中，通常會(huì)繞著自身體軸高速旋轉(zhuǎn)來(lái)提高飛行穩(wěn)定性，增加飛行軌跡的精度。這種錐進(jìn)運(yùn)動(dòng)在保證了飛行性能的同時(shí)也帶來(lái)了負(fù)面效應(yīng)，當(dāng)彈體存在攻角時(shí)，由于旋轉(zhuǎn)引起的彈體表面流場(chǎng)相對(duì)攻角平面不對(duì)稱，產(chǎn)生垂直于該平面的力和繞重心的力矩。該現(xiàn)象由德國(guó)科學(xué)家馬格努斯首先發(fā)現(xiàn)，故將該力和力矩稱為馬格努斯力和力矩[1-3]。馬格努斯效應(yīng)產(chǎn)生的主要原因有二：一是旋轉(zhuǎn)引起的附面層畸變，即附面層位移厚度沿著彈體周向發(fā)生變化；二是離心力因素。由于馬格努斯效應(yīng)影響很大，因此準(zhǔn)確預(yù)測(cè)其力和力矩，減小其對(duì)飛行的影響，是當(dāng)前研究彈體流場(chǎng)特性的熱點(diǎn)之一。

目前，使用CFD方法開(kāi)展彈體馬格努斯力矩效應(yīng)數(shù)值計(jì)算已經(jīng)成為一種常用方法，國(guó)內(nèi)外也取得了很多進(jìn)展和成果。國(guó)外，W.B.Sturek等[4]基于薄層假設(shè)，對(duì)尖錐等旋成體的馬格努斯效應(yīng)進(jìn)行了數(shù)值計(jì)算分析，研究了幾何參數(shù)的變化對(duì)彈丸氣動(dòng)特性的影響；C.J.Nietubicz等[5]基于動(dòng)態(tài)網(wǎng)格技術(shù)，對(duì)跨聲速?gòu)椡璧鸟R格努斯效應(yīng)進(jìn)行了數(shù)值模擬，分析了不同參數(shù)對(duì)該效應(yīng)的影響；M.Muto等[6]基于數(shù)值方法對(duì)反馬格努斯效應(yīng)引起的邊界層轉(zhuǎn)捩現(xiàn)象進(jìn)行研究，揭示了高速旋轉(zhuǎn)的彈丸的邊界層流動(dòng)細(xì)節(jié)；D.Klatt等[7]則對(duì)超聲速飛行的旋轉(zhuǎn)彈丸的馬格努斯現(xiàn)象進(jìn)行了研究，使用的方法為非結(jié)構(gòu)網(wǎng)格動(dòng)態(tài)重構(gòu)模擬旋轉(zhuǎn)運(yùn)動(dòng)，對(duì)比試驗(yàn)結(jié)果分析了激波現(xiàn)象；S.I.Silton[8]對(duì)彈丸的亞跨超聲速時(shí)的馬格努斯效應(yīng)進(jìn)行了模擬，馬格努斯力及力矩結(jié)果相比升阻力系數(shù)較試驗(yàn)值偏差較大，說(shuō)明旋轉(zhuǎn)效應(yīng)的模擬精度還有待提升；自2004年開(kāi)始，J.Despirito[9]使用DES模型開(kāi)展了M910旋轉(zhuǎn)彈丸的馬格努斯力矩效應(yīng)研究，得到了更為精確的計(jì)算結(jié)果。國(guó)內(nèi)，盡管起步較晚，但進(jìn)展迅速。陳東陽(yáng)等[10]利用Fluent軟件的動(dòng)網(wǎng)格技術(shù)，對(duì)帶有結(jié)構(gòu)誤差的M910彈丸的馬格努斯效應(yīng)進(jìn)行了數(shù)值計(jì)算，分析了質(zhì)量不對(duì)稱引起的馬格努斯力偏差；劉周等[11]利用DES方法對(duì)彈丸在大攻角時(shí)的馬格努斯效應(yīng)進(jìn)行了計(jì)算分析；雷娟棉等[12]基于滑移網(wǎng)格方法、馬杰等[13]使用滑移動(dòng)網(wǎng)格方法對(duì)SOCBT彈丸的氣動(dòng)特性進(jìn)行了研究，提出并驗(yàn)證了適用于小攻角、超聲速情形的彈丸馬格努斯力和力矩系數(shù)經(jīng)驗(yàn)公式。

上述研究均基于較復(fù)雜的計(jì)算方法，例如動(dòng)網(wǎng)格、滑移網(wǎng)格等，盡管精度較高，但分析效率相對(duì)較低，限制了這些方法在彈體快速設(shè)計(jì)中的應(yīng)用。

本文針對(duì)該問(wèn)題，使用一種簡(jiǎn)化的定常求解方法——旋轉(zhuǎn)壁面法，進(jìn)行彈體馬格努斯效應(yīng)的快速模擬，并與試驗(yàn)值進(jìn)行對(duì)比分析；針對(duì)影響馬格努斯效應(yīng)的附面層堆積，開(kāi)展附面層網(wǎng)格節(jié)點(diǎn)數(shù)和不同計(jì)算模型的適應(yīng)性研究，對(duì)比分析這些因素對(duì)本文方法的影響。

1 計(jì)算方法

1.1 常用方法

為了計(jì)算彈體的旋轉(zhuǎn)效應(yīng)，常用的方法有三種：①基于剛性動(dòng)網(wǎng)格技術(shù)的非定常計(jì)算方法，②基于滑移網(wǎng)格的非定常計(jì)算方法，③基于滑移網(wǎng)格的準(zhǔn)定常計(jì)算方法。其中，基于剛性動(dòng)網(wǎng)格的非定常計(jì)算方法將整個(gè)流場(chǎng)域視為剛性的整體，網(wǎng)格跟隨彈體的旋轉(zhuǎn)而轉(zhuǎn)動(dòng)，這樣處理的好處是網(wǎng)格的拓?fù)浣Y(jié)構(gòu)不因旋轉(zhuǎn)而發(fā)生改變，但是外邊界處會(huì)附加旋轉(zhuǎn)速度，該方法如圖1所示；基于滑移網(wǎng)格的非定常方法將流場(chǎng)域劃分為動(dòng)、靜兩部分，動(dòng)域按照設(shè)定的旋轉(zhuǎn)方式進(jìn)行運(yùn)動(dòng)，動(dòng)靜域之間通過(guò)交界面進(jìn)行數(shù)據(jù)傳遞，動(dòng)域的網(wǎng)格整體運(yùn)動(dòng)，因此計(jì)算網(wǎng)格質(zhì)量較高，該方法也是一種較為精確的計(jì)算方法，如圖2所示；基于滑移網(wǎng)格的準(zhǔn)定常計(jì)算方法是將旋轉(zhuǎn)參考系附加到滑移網(wǎng)格上，此時(shí)動(dòng)域的網(wǎng)格并不會(huì)真正旋轉(zhuǎn)，這種準(zhǔn)定常方法是一種簡(jiǎn)化的高效方法，如圖3所示。

圖1 剛性動(dòng)網(wǎng)格

圖2 滑移動(dòng)網(wǎng)格

圖3 旋轉(zhuǎn)參考系方法

1.2 基于旋轉(zhuǎn)壁面法的馬格努斯效應(yīng)計(jì)算方法

圖4 旋轉(zhuǎn)壁面法

本文計(jì)算全部使用ANSYS FLUENT軟件完成。為了提高計(jì)算速度和收斂性，計(jì)算時(shí)先進(jìn)行無(wú)旋工況的模擬，此時(shí)的邊界條件設(shè)置如圖5所示。在此基礎(chǔ)上設(shè)定彈體壁面的旋轉(zhuǎn)速度和轉(zhuǎn)軸，繼續(xù)計(jì)算得到此工況的馬格努斯力和力矩，旋轉(zhuǎn)情形的邊界條件設(shè)置如圖6所示。

圖5 無(wú)旋轉(zhuǎn)情形邊界條件

圖6 旋轉(zhuǎn)情形邊界條件

旋轉(zhuǎn)壁面法并不是一種很新的方法，在計(jì)算流體領(lǐng)域多有應(yīng)用，但是在多對(duì)稱面彈體的馬格努斯力矩?cái)?shù)值計(jì)算方面卻鮮少涉及。本文將該方法引入旋轉(zhuǎn)彈體的流場(chǎng)計(jì)算中，形成一種簡(jiǎn)化的快速馬格努斯力矩計(jì)算方法。需要注意的是，該快速計(jì)算方法主要用來(lái)計(jì)算多對(duì)稱面的旋成體馬格努斯力矩效應(yīng)，并不適用于不規(guī)則的外形。

1.3 數(shù)值方法

采用的控制方程為三維積分形式的雷諾平均N-S方程：

(1)

式中：V為任意控制體；W為守恒變量；F為無(wú)粘(對(duì)流)通矢量項(xiàng)；Fv為粘性通量；?V為控制體的邊界；n為控制體邊界單位外法向矢量；Re為計(jì)算的雷諾數(shù)。

空間離散采用二階迎風(fēng)格式——通量差分分裂(Roe-FDS)格式，采用隱式LU-SGS(Lower-Upper Symmetric Gauss-Seidel)格式迭代求解，應(yīng)用當(dāng)?shù)貢r(shí)間步長(zhǎng)、殘值光順、預(yù)處理和多重網(wǎng)格加速收斂。湍流模型采用對(duì)逆壓梯度流動(dòng)模擬精度較高的k-ωSST(剪切應(yīng)力輸運(yùn))湍流模型。

2 計(jì)算模型及網(wǎng)格

2.1 計(jì)算模型

計(jì)算模型為一6倍口徑彈丸，其名稱是Secant-Ogive-Clinder with Boat Tail，簡(jiǎn)稱SOCBT[14]。該模型是評(píng)估馬格努斯力矩效應(yīng)的標(biāo)準(zhǔn)模型，由正割風(fēng)帽頭部、圓柱彈身、船尾形尾部三部分構(gòu)成，船尾收縮角β=7°。其幾何模型和尺寸如圖7所示。

圖7 SOCBT彈丸模型

計(jì)算坐標(biāo)系如圖8所示。θ=0°為迎風(fēng)面，θ=180°為背風(fēng)面。z向的力Cz和繞著y軸的力矩Cmy即為馬格努斯力和力矩，其參考點(diǎn)位于彈體頭部原點(diǎn)處，參考面積為彈體的最大橫截面積Sref=0.002 57m2，參考長(zhǎng)度為彈體的最大直徑Cref=d=0.057 20 m。

在以往的發(fā)展中，由于旅游業(yè)在一些項(xiàng)目中利潤(rùn)較高，從業(yè)者人數(shù)大量的增加，在一定的程度上使得旅游業(yè)內(nèi)部的環(huán)境受到了沖擊，人們普遍的對(duì)于旅游公司以及導(dǎo)游等具有不佳的印象。在大數(shù)據(jù)時(shí)代到來(lái)之后，人們?cè)谏畹倪^(guò)程中，能夠以更加透明的方式對(duì)于旅游業(yè)內(nèi)部的情況進(jìn)行了解。并且在此基礎(chǔ)上，無(wú)形對(duì)于行業(yè)的內(nèi)部構(gòu)成了一種監(jiān)督的機(jī)制，使得行業(yè)內(nèi)部的環(huán)境得到了凈化，服務(wù)的質(zhì)量有所提升。

圖8 使用的坐標(biāo)系

2.2 計(jì)算網(wǎng)格

計(jì)算網(wǎng)格使用ANSYS ICEM CFD軟件生成的高質(zhì)量結(jié)構(gòu)網(wǎng)格，其拓?fù)浣Y(jié)構(gòu)如圖9所示，表面網(wǎng)格和截面網(wǎng)格分別如圖10～圖11所示?？偩W(wǎng)格量540萬(wàn)，附面層首層網(wǎng)格高度y+≤1.0。

圖9 網(wǎng)格拓?fù)浣Y(jié)構(gòu)

圖10 整體計(jì)算網(wǎng)格

圖11 壁面網(wǎng)格

3 算例與驗(yàn)證分析

3.1 無(wú)旋轉(zhuǎn)驗(yàn)證

首先針對(duì)無(wú)旋工況進(jìn)行計(jì)算，來(lái)驗(yàn)證本文數(shù)值方法的有效性。計(jì)算狀態(tài)為：來(lái)流馬赫數(shù)Ma=3，大氣溫度T=310 K，大氣壓強(qiáng)P=298 500 Pa，攻角α=6.34°。計(jì)算得到的彈體對(duì)稱面以及截面的壓力云圖分別如圖12～圖13所示。

圖12 對(duì)稱面壓力云圖

圖13 截面x/d=5.77處的壓力云圖

從圖12可以看出：超聲速的氣流在彈體周?chē)纬啥嗟兰げǎ⒃诖矃^(qū)域形成膨脹波，物理上符合流場(chǎng)的基本規(guī)律。

從圖13可以看出：無(wú)旋時(shí)的截面壓力分布完全對(duì)稱。

x/d=4.88和x/d=5.77截面處的詳細(xì)壓力分布分別如圖14～圖15所示。

圖14 x/d=4.88處的壓力對(duì)比

圖15 x/d=5.77處的壓力對(duì)比

從圖14～圖15可以看出：兩個(gè)截面上的壓力計(jì)算值與試驗(yàn)值吻合得很好，表明本文的數(shù)值方法是有效的。

3.2 旋轉(zhuǎn)時(shí)的馬格努斯效應(yīng)

在無(wú)旋驗(yàn)證的基礎(chǔ)上，使用所建立的旋轉(zhuǎn)壁面法進(jìn)行彈體的馬格努斯力矩效應(yīng)的數(shù)值模擬。計(jì)算的攻角為0°、3°、5°、6.34°、7°、10°、12°以及15°。為了保證計(jì)算條件與動(dòng)態(tài)試驗(yàn)一致，取無(wú)量綱的旋轉(zhuǎn)角速度Ω*=Ωd/V=0.19，其中Ω為實(shí)際的旋轉(zhuǎn)角速度。計(jì)算得到的各個(gè)攻角下的馬格努斯力和力矩系數(shù)分別如圖16～圖17所示。

圖16 馬格努斯力系數(shù)

圖17 馬格努斯力矩系數(shù)

從圖16～圖17可以看出：在中小攻角范圍，計(jì)算值與試驗(yàn)值吻合較好，隨著攻角的增大，彈體背風(fēng)區(qū)出現(xiàn)分離，而本文使用的湍流模型對(duì)分離的捕捉能力有限，此時(shí)計(jì)算結(jié)果與試驗(yàn)值的偏差變大。但總體來(lái)看，本文的計(jì)算方法是有效的。

3.3 數(shù)值計(jì)算方法影響因素分析

3.3.1 附面層的影響

馬格努斯效應(yīng)產(chǎn)生的主要原因之一是旋轉(zhuǎn)引起的附面層堆積，因此對(duì)周向的附面層差異進(jìn)行足夠精確地模擬是獲取準(zhǔn)確的馬格努斯力和力矩的前提。為了準(zhǔn)確地描述附面層，需要在附面層內(nèi)部布置足夠的網(wǎng)格節(jié)點(diǎn)。分別設(shè)定附面層的網(wǎng)格節(jié)點(diǎn)數(shù)量為21、61和101，對(duì)旋轉(zhuǎn)時(shí)的流場(chǎng)進(jìn)行計(jì)算。使用的湍流模型仍為k-ωSST模型，以初始攻角α=6.34°的工況為例進(jìn)行分析。三種模型計(jì)算得到的馬格努斯力矩系數(shù)如表1所示。

表1 不同附面層節(jié)點(diǎn)數(shù)計(jì)算結(jié)果對(duì)比

從表1可以看出：隨著附面層節(jié)點(diǎn)數(shù)的增加，計(jì)算結(jié)果與試驗(yàn)值吻合得越好，表明更多的附面層節(jié)點(diǎn)能夠?qū)Χ逊e效應(yīng)做出更準(zhǔn)確地模擬；而網(wǎng)格節(jié)點(diǎn)為61時(shí)，其計(jì)算結(jié)果與節(jié)點(diǎn)為101的計(jì)算結(jié)果差別較小，從減小計(jì)算代價(jià)的角度考慮，可以選擇節(jié)點(diǎn)數(shù)為61進(jìn)行模擬。

3.3.2 湍流模型的影響

旋轉(zhuǎn)引起的馬格努斯力矩不僅導(dǎo)致了附面層的堆積，而且由于粘性致使堆積區(qū)域的分離也直接影響著馬格努斯力的準(zhǔn)確計(jì)算。尤其是隨著初始攻角的增大，背風(fēng)區(qū)的分離區(qū)域擴(kuò)大，與堆積區(qū)的分離相互干擾，導(dǎo)致整個(gè)彈體的流場(chǎng)表現(xiàn)為較復(fù)雜的狀態(tài)。而不同的湍流模型對(duì)這種復(fù)雜流態(tài)的捕捉能力不同，本文選擇一方程的S-A模型、兩方程的k-ωSST模型和SAS自適應(yīng)尺度模型[15]進(jìn)行6.34°攻角和15°攻角的對(duì)比計(jì)算。其中SAS模型使用類似于DES模型的尺度因子，能夠?qū)Ψ蛛x現(xiàn)象做出較為精確的模擬。

6.34°攻角時(shí)的計(jì)算結(jié)果如表2所示，可以看出：由于此時(shí)的分離并不特別嚴(yán)重，三種模型的計(jì)算結(jié)果差別不明顯。而15°攻角時(shí)的馬格努斯力和力矩的對(duì)比如表3所示，可以看出：此時(shí)差別比較明顯，流場(chǎng)背風(fēng)區(qū)出現(xiàn)大范圍的分離流動(dòng)，且與激波形成復(fù)雜的干擾；使用SAS模型得到的結(jié)果顯然更接近試驗(yàn)值，主要是因?yàn)樵撃Ｐ拖啾纫话愕耐牧髂Ｐ湍軌虿蹲降礁鼮榫?xì)、更為準(zhǔn)確的分離效應(yīng)。

表2 6.34°攻角下不同計(jì)算模型的結(jié)果

表3 15°攻角下不同計(jì)算模型的結(jié)果

不同湍流模型的計(jì)算表明，為了能夠更好地獲取彈體的馬格努斯力矩效應(yīng)，對(duì)于不同攻角的模擬需要選擇不同的計(jì)算模型，達(dá)到兼顧效率和精度的效果。

4 結(jié) 論

(1) 本文提出的基于旋轉(zhuǎn)壁面法的多對(duì)稱面彈體旋轉(zhuǎn)效應(yīng)數(shù)值模擬方法，不使用動(dòng)網(wǎng)格和多參考系方法，將旋轉(zhuǎn)速度附加于壁面上，使用定常方法進(jìn)行計(jì)算，效率較高。對(duì)旋轉(zhuǎn)彈體的馬格努斯力和力矩系數(shù)的計(jì)算結(jié)果驗(yàn)證了方法的有效性。

(2) 馬格努斯力矩效應(yīng)的起因之一是旋轉(zhuǎn)引起的附面層堆積，為了準(zhǔn)確地捕捉這種效應(yīng)，需要在附面層區(qū)域布置足夠的網(wǎng)格節(jié)點(diǎn)。

(3) 大攻角時(shí)，旋轉(zhuǎn)引起的附面層堆積同分離區(qū)域相互干擾，對(duì)于超聲速的彈體來(lái)講，還有多重激波的相互摻雜，這種復(fù)雜的非線性流場(chǎng)使用一般的湍流模型很難精確模擬，需要借助更高精度的模型，SAS模型作為一種新的大攻角計(jì)算方法，使用了類似DES方法的尺度因子，能夠在降低網(wǎng)格要求的情形下得到與DES類方法接近的計(jì)算效果。

[1] Seifert J. A review of the Magnus effect in aeronautics[J]. Progress in Aerospace Sciences, 2012, 55(5): 17-45.

[2] Zhao J, Hou Q, Jin H, et al. Discussion on improving Magnus effect of cylinder based on CFD[C]∥IEEE International Conference on Mechatronics and Automation. IEEE, 2013: 539-543.

[3] Cayzac R, Carette E, Denis P, et al. Magnus effect: physical origins and numerical prediction[J]. Journal of Applied Mechanics, 2011, 78(5): 748-760.

[4] Sturek W B, Kayser L D, Nietubicz C J, et al. Computations of Magnus effects for a yawed, spinning body of revolution[J]. AIAA Journal, 2012, 16(7): 687-692.

[5] Nietubicz C J, Sturek W B, Heavey K R. Computations of projectile Magnus effect at transonic velocities[J]. AIAA Journal, 2012, 23(7): 15.

[6] Muto M, Tsubokura M, Oshima N. Numerical visualization of boundary layer transition when negative Magnus effect occurs[J]. Journal of Visualization, 2012, 15(3): 261-268.

[7] Klatt D, Hruschka R, Leopold F. Numerical and experimental investigation of the Magnus effect in supersonic flows[C]∥AIAA Applied Aerodynamics Conference, 2013.

[8] Silton S I. Navier-Stokes computations for a spinning projectile from subsonic to supersonic speeds[J]. Journal of Spacecraft and Rockets, 2005, 42(2): 223-231.

[9] Despirito J. CFD prediction of Magnus effects in subsonic to supersonic flow[R]. AIAA-2008-0427, 2008.

[10] 陳東陽(yáng), Laith K A, 芮筱亭, 等. 結(jié)構(gòu)誤差對(duì)旋轉(zhuǎn)穩(wěn)定彈丸氣動(dòng)特性影響的數(shù)值模擬[J]. 空氣動(dòng)力學(xué)學(xué)報(bào), 2014, 32(5): 705-711.

Chen Dongyang, Laith K A, Rui Xiaoting, et al. Numerical simulation of a spinning stabilized projectile aerodynamic characteristics effected by structure errors[J]. Acta Aerodynamica Sinica, 2014, 32(5): 705-711.(in Chinese)

[11] 劉周, 謝立軍, 楊云軍, 等. 彈丸旋轉(zhuǎn)空氣動(dòng)力效應(yīng)非定常數(shù)值模擬[J]. 航空學(xué)報(bào), 2016, 37(5): 1401-1410.

Liu Zhou, Xie Lijun, Yang Yunjun, et al. Unsteady numerical simulation of aerodynamic effect of a spinning projectile[J]. Acta Aeronautica et Astronautica Sinica, 2016, 37(5): 1401-1410.(in Chinese)

[12] 雷娟棉, 李田田, 黃燦. 高速旋轉(zhuǎn)彈丸馬格努斯效應(yīng)數(shù)值研究[J]. 兵工學(xué)報(bào), 2013, 34(6): 718-725.

Lei Juanmian, Li Tiantian, Huang Can. A numerical investigation of Magnus effect for high-speed spinning projectile[J]. Acta Armamenyarii, 2013, 34(6): 718-725.(in Chinese)

[13] 馬杰, 陳志華, 姜孝海. 高速旋轉(zhuǎn)條件下的彈丸氣動(dòng)特性研究[J]. 彈道學(xué)報(bào), 2015, 27(2): 1-6.

Ma Jie, Chen Zhihua, Jiang Xiahai. Aerodynamic characteristics of a projectile with high spinning speeds[J]. Journal of Ballistics, 2015, 27(2): 1-6.(in Chinese)

[14] Onn S C, Su A, Wei C K, et al. Computational drag and Magnus force reduction for a transonic spinning projectile using passive porosity[J]. Computer Methods on Applied Mechanics and Engineering, 2001, 190(46/47): 6125-6139.

[15] Menter F, Egorov Y. The scale-adaptive simulation method for unsteady turbulent flow predictions. Part 1: theory and model description[J]. Flow Turbulence and Combustion, 2010, 85(1): 113-138.