張 超，李 丹，田 磊，楊文鳳，王學(xué)德，馬庚軍

(1.中國飛行試驗(yàn)研究院, 西安 710089; 2.南京理工大學(xué) 能源與動(dòng)力工程學(xué)院, 南京 210094)

為了大幅提高不帶控制系統(tǒng)彈箭的命中精度和穩(wěn)定性，一般都采用彈箭低速旋轉(zhuǎn)的方法[1]，而實(shí)現(xiàn)彈箭旋轉(zhuǎn)的方式之一是采用傾斜尾翼，但尾翼傾斜角選取得不合理，輕則導(dǎo)致彈箭飛行不穩(wěn)定，重則導(dǎo)致掉彈危險(xiǎn)。基于此背景，研究尾翼傾斜角度的大小對其氣動(dòng)特性影響就變得非常有意義。

DeSpirito[2]采用小擾動(dòng)法對彈丸的Magnus效應(yīng)進(jìn)行了計(jì)算，求出Magnus力和力矩。Klatt等[3]對比風(fēng)洞試驗(yàn)結(jié)果和采用RANS方法下的計(jì)算結(jié)果，得到了大攻角下Magnus效應(yīng)的變化規(guī)律。Leroy[4]對舵身組合彈箭的氣動(dòng)特性計(jì)算分析表明，在大攻角下，其滾轉(zhuǎn)特性和Magnus特性參數(shù)呈非線性變化趨勢。Nietubicz[5]基于薄層假設(shè)的非定常Navier-Stokes方程發(fā)現(xiàn)邊界層的非對稱畸變對彈丸Magnus力的貢獻(xiàn)最大。

2005年，高旭東等[6]比較了穩(wěn)態(tài)和瞬態(tài)方法對旋轉(zhuǎn)彈流場的計(jì)算結(jié)果，發(fā)現(xiàn)兩者計(jì)算得到的動(dòng)導(dǎo)數(shù)相差較小，但是瞬態(tài)計(jì)算時(shí)間遠(yuǎn)大于穩(wěn)態(tài)時(shí)間。2011年，周德娟[7]利用旋轉(zhuǎn)坐標(biāo)法模擬了不同尾翼數(shù)下旋轉(zhuǎn)彈的滾轉(zhuǎn)特性，發(fā)現(xiàn)尾翼數(shù)目與滾轉(zhuǎn)力矩系數(shù)呈反比關(guān)系。張超等[8]計(jì)算了不同轉(zhuǎn)速下彈箭的氣動(dòng)特性，發(fā)現(xiàn)Magnus力隨滾轉(zhuǎn)角呈正弦變化趨勢，轉(zhuǎn)速與滾轉(zhuǎn)力矩呈線性變化關(guān)系。石磊等[9]基于修正一方程Spalart-Allmaras湍流模型和兩方程SST湍流模型計(jì)算了高速旋轉(zhuǎn)的ANF(Army-Navy Finner)標(biāo)模的氣動(dòng)力及力矩。

進(jìn)入21世紀(jì)以來，各國科研工作者對旋轉(zhuǎn)彈展開了系統(tǒng)而深入的研究，已經(jīng)掌握Magnus效應(yīng)、滾轉(zhuǎn)特性等氣動(dòng)參數(shù)的變化機(jī)理，但對實(shí)現(xiàn)彈箭旋轉(zhuǎn)的方式——尾翼傾斜角對彈箭氣動(dòng)參數(shù)的影響，這方面的研究文獻(xiàn)非常少。基于此背景，本文提出基于滑移網(wǎng)格技術(shù)的旋轉(zhuǎn)彈氣動(dòng)特性數(shù)值仿真方法，系統(tǒng)地研究尾翼傾斜角對彈箭氣動(dòng)特性的影響。

1 數(shù)值模擬方法

滑移網(wǎng)格技術(shù)具有可信度高，運(yùn)行速度快等優(yōu)勢[9]，普遍地應(yīng)用在彈箭旋轉(zhuǎn)的氣動(dòng)計(jì)算方面，圖1為滑移網(wǎng)格區(qū)域劃分示意圖。

圖1 滑移網(wǎng)格區(qū)域劃分示意圖

控制方程采用非定常守恒型Navier-Stokes方程組

可實(shí)現(xiàn)k-ε湍流模型能較為準(zhǔn)確地捕捉到旋轉(zhuǎn)流動(dòng)的流場細(xì)節(jié)[10-11]，因此本文采用該模型進(jìn)行數(shù)值仿真，空間離散采用有限體積法，對流項(xiàng)設(shè)置為AUSM+格式，黏性項(xiàng)設(shè)置為中心差分格式。

對于非定常計(jì)算結(jié)果影響較大關(guān)鍵參數(shù)之一是時(shí)間步長，合理的時(shí)間步長可以加速收斂所需的時(shí)間，因此本文采用雙時(shí)間步長技術(shù)。

2 算法的有效性驗(yàn)證

2.1 驗(yàn)證模型

為了驗(yàn)證該數(shù)值仿真算法的有效性，本文選取了F4彈箭模型[7]，該模型的具體尺寸如圖2所示，其計(jì)算工況有關(guān)數(shù)據(jù)如表1所示，無量綱轉(zhuǎn)速的計(jì)算公式為彈徑(轉(zhuǎn)速/聲速，N為迭代次數(shù)，Δt為計(jì)算步長。

表1 計(jì)算工況

升力系數(shù)、阻力系數(shù)、繞過某一參考軸的力矩系數(shù)分別為

圖2 F4彈箭模型

2.2 算法的有效性驗(yàn)證

圖3給出了F4模型的氣動(dòng)系數(shù)隨Ma值的變化曲線，并給出了相同工況下參考文獻(xiàn)的計(jì)算結(jié)果。由圖可知：文本中的計(jì)算結(jié)果與參考文獻(xiàn)中計(jì)算值高度吻合，誤差均在5%以內(nèi)，因此可以說明本文選取的數(shù)值計(jì)算方法適用于全Ma數(shù)條件下，并且有較高的精度。

圖3 氣動(dòng)系數(shù)隨Ma值的變化曲線

3 彈箭Magnus效應(yīng)數(shù)值計(jì)算

3.1 模型、網(wǎng)格和計(jì)算狀態(tài)

計(jì)算模型如圖4(a)，其尺寸參數(shù)如圖4(b)所示。彈徑Dm=412 mm，重心Xmc=4.2D，彈長Dl=7.2D，翼根弦長br=1.1D，翼梢弦長bt=0.4D，網(wǎng)格模型如圖5。選取的計(jì)算狀態(tài)參數(shù)如表2所示。

圖4 模型及尺寸

圖5 彈箭網(wǎng)格模型

表2 計(jì)算狀態(tài)參數(shù)

3.2 網(wǎng)格和時(shí)間步長無關(guān)性驗(yàn)證

為了選取合適的時(shí)間步長和網(wǎng)格數(shù)，本文給出了Ma=1.6、α=2°、γ=1°工況下，時(shí)間步長分別為0.004、0.006、0.008以及0.01下各氣動(dòng)系數(shù)，如圖 6所示。可見隨著時(shí)間步長的不斷減小，其氣動(dòng)特性系數(shù)趨于收斂，因此選取后續(xù)計(jì)算時(shí)間步長為0.006。

圖6 不同時(shí)間步長下的氣動(dòng)特性系數(shù)曲線

在網(wǎng)格數(shù)M分別等于8×106、7×106、6×106、5×106下氣動(dòng)系數(shù)如表 3所示?？梢婋S著網(wǎng)格數(shù)的不斷增加，相鄰兩組網(wǎng)格數(shù)下的氣動(dòng)力或力矩的差值越來越小，表明該系數(shù)逐漸穩(wěn)定。由于最后兩組網(wǎng)格數(shù)下的計(jì)算結(jié)果差異較小(四組氣動(dòng)特性系數(shù)的差值分別為0.11%、0.51%、0.72%、0.76%)，所以后續(xù)計(jì)算選取網(wǎng)格數(shù)7×106。

表3 網(wǎng)格數(shù)對氣動(dòng)特性系數(shù)

3.3 不同狀態(tài)下的流場結(jié)構(gòu)

圖7給出了在Ma=2.0、α=2°時(shí)尾翼處的壓力云圖，其中x表示截面位置，L表示彈長，圖中的1～4分別代表4片尾翼。從圖 7(a)可知：彈體附近的壓力較其他部分的壓力高，并且呈“蝴蝶”形狀分布；由于此時(shí)彈箭關(guān)于xy平面不對稱，并且由于攻角產(chǎn)生的橫流的影響，導(dǎo)致尾翼兩側(cè)的壓力分布不對稱。同時(shí)發(fā)現(xiàn)尾翼4兩側(cè)的壓力差最大。這是由于尾翼3并未完全擋住尾翼4，尾翼4的下翼面仍然處于橫流(即迎風(fēng)面)中，上翼面完全處于背風(fēng)區(qū)。尾翼1兩側(cè)的壓力差最小，這是由于尾翼2完全遮擋了尾翼1，導(dǎo)致尾翼1完全處于背風(fēng)區(qū)。對比圖 7(a)、(b)、(c)發(fā)現(xiàn)：隨著尾翼傾斜角的增大，尾翼1處兩側(cè)的壓力差逐漸增大，這是由于隨著尾翼傾斜角的增大，尾翼兩側(cè)的氣流分布不均進(jìn)一步加劇。

圖7 不同尾翼傾斜角下尾翼處(x/L=0.95)壓力云圖(Ω=0.003 8)

圖8給出了在Ma=2.0、α=2°時(shí)彈體不同截面的渦量強(qiáng)度，其中x表示截面位置，L表示彈長。同時(shí)為了更好地顯示渦量強(qiáng)度，將渦量強(qiáng)度大于10萬以及小于1 500的區(qū)域不予顯示。由圖 8(a)可知：盡管尾翼傾斜角為0，但是由于轉(zhuǎn)速不為0且此時(shí)彈箭關(guān)于攻角平面不對稱，因此在彈體表面以及尾翼附近產(chǎn)生渦流。對比圖 8(a)、(b)、(c)可知：隨著尾翼傾斜角的增大，氣流由于混合產(chǎn)生的擾動(dòng)加劇，從而使得渦流強(qiáng)度加劇。

圖8 不同尾翼傾斜角下尾翼處(x/L=0.95)的渦量強(qiáng)度圖(Ω=0.003 8)

3.4 不同狀態(tài)下氣動(dòng)特性分析

本文中彈箭模型的尾翼采用差動(dòng)布置方式，即對一對尾翼而言，一片尾翼的傾斜角為+γ，則另一片尾翼的傾斜角為-γ。

從圖 9可知，尾翼傾斜角從0°增加到2°的過程中，平均阻力系數(shù)處于非定比值增大趨勢。這主要是因?yàn)閷τ谝粚Σ顒?dòng)布置的尾翼而言，一側(cè)尾翼的局部附加攻角增加而另一側(cè)尾翼的局部附加攻角減小，并且阻力系數(shù)隨攻角的增大呈現(xiàn)非線性增加特性(近似拋物線增長規(guī)律)，即由+γ引起的阻力增量遠(yuǎn)大于由-γ引起的減小量，正尾翼傾斜角引起的阻力增加量減去負(fù)尾翼傾斜角引起的阻力減小量所得的結(jié)果即為整個(gè)彈箭的阻力增加量，并且γ越大則平均阻力系數(shù)增加越明顯，由表4可知，在同一Ma數(shù)下，彈體提供的阻力基本保持不變，而尾翼提供的阻力隨尾翼傾斜角的增大而顯著增大；表5給出了尾翼傾斜角γ引起的平均阻力系數(shù)增量。

在小攻角范圍內(nèi)，升力系數(shù)呈定比值(即線性變化)變化規(guī)律，因此由+γ，-γ引起的升力增大量和減小量大小相等方向相反，因此整個(gè)彈箭總升力系數(shù)隨著尾翼傾斜角的變化基本保持不變。

圖9 平均阻力系數(shù)和升力系數(shù)隨Ma數(shù)的變化曲線

表4 4°攻角下不同部件對彈箭阻力增加量

表5 尾翼傾斜角γ引起的平均阻力系數(shù)增加量 %

從圖 10(a)可知，隨著尾翼傾斜角從0°逐漸增大到2°，三條平均俯仰力矩系數(shù)曲線基本重合，產(chǎn)生該現(xiàn)象的原因是在小攻角情況下，升力在彈箭上的作用點(diǎn)不隨攻角變化而變化，因此升力減小量ΔL左對質(zhì)心取矩形成的俯仰力矩ΔML-左和升力增加量ΔL右對質(zhì)心形成的俯仰力矩ΔML-右大小相等且方向相反，總俯仰力矩系數(shù)隨著尾翼傾斜角的變化基本保持不變。

從圖 10(b)可知，隨著尾翼傾斜角從0°增加到2°，平均滾轉(zhuǎn)力矩系數(shù)呈增大趨勢，且當(dāng)尾翼傾斜角為0時(shí)，該系數(shù)為負(fù)值，即為滾轉(zhuǎn)阻尼力矩系數(shù)，這是由于氣流黏性的存在，阻礙彈箭的旋轉(zhuǎn)。當(dāng)傾斜角大于0°時(shí)，該系數(shù)均為正值，這是由于空氣動(dòng)力引起的滾轉(zhuǎn)力矩始終大于空氣黏性引起的滾轉(zhuǎn)阻尼力矩的緣故。由表 6可知，在同一攻角和Ma數(shù)下，隨著尾翼傾斜角的增加，平均滾轉(zhuǎn)力矩系數(shù)典型的定比值呈增大趨勢，這是由于彈體左側(cè)翼面產(chǎn)生一個(gè)向下的附加升力ΔL左而彈體右側(cè)翼面產(chǎn)生一個(gè)向上的附加升力ΔL右，這兩個(gè)附加升力對OX軸取矩則形成了滾轉(zhuǎn)力矩，并且由于該力矩方向與參考方向相同，所以平均滾轉(zhuǎn)力矩系數(shù)為正值。在小攻角時(shí)，升力隨攻角呈線性變化規(guī)律，因此由附加升力產(chǎn)生的滾轉(zhuǎn)力矩隨γ的增大呈線性變化規(guī)律。

圖10 平均俯仰力矩和滾轉(zhuǎn)力矩系數(shù)隨Ma數(shù)的變化曲線

表6 尾翼傾斜角引起的滾轉(zhuǎn)力矩系數(shù)增加量

由圖 11可知，平均Magnus力系數(shù)隨Mach數(shù)呈現(xiàn)先增大后下降的趨勢，并且在Ma=1.1附近取得最大值；平均Magnus力系數(shù)隨著尾翼傾斜角的變化呈現(xiàn)非定比值增大趨勢；平均Magnus力矩系數(shù)的變化趨勢與該系數(shù)的變化趨勢一致，不再贅述。

由表7可知：對比兩者所占百分比可知：全彈的Magnus力系數(shù)主要由壓力分量貢獻(xiàn)，尤其是Ma=1.1時(shí)；黏性分量對全彈的平均Magnus力系數(shù)的貢獻(xiàn)可以忽略不計(jì)。

4 結(jié)論

1) 隨著尾翼傾斜角的增加，平均阻力系數(shù)呈現(xiàn)非定比值增大趨勢，平均滾轉(zhuǎn)力矩系數(shù)呈定比值增大趨勢；

2) 隨著尾翼傾斜角的增加，平均升力及俯仰力矩系數(shù)不變。

3) 隨著尾翼傾斜角的增大，平均Magnus力和力矩系數(shù)呈非線性增大趨勢，并且相比于壓差分量，黏性分量對Magnus力的貢獻(xiàn)可忽略不計(jì)，尤其在臨界Ma數(shù)附近。