滾轉(zhuǎn)和錐進運動對彈箭動導(dǎo)數(shù)求解的影響

陳 亮, 劉榮忠,*, 郭 銳, 陳福紅, 侯俊亮, 楊永亮, 邢柏陽, 高 科

(1. 南京理工大學(xué) 機械工程學(xué)院, 江蘇 南京 290014; 2. 中國航天科技集團公司第七研究院 第七設(shè)計部, 四川 成都 610100)

0 引 言

動穩(wěn)定性導(dǎo)數(shù)是影響彈箭飛行穩(wěn)定性和可控性的重要氣動參數(shù)，因此在彈箭氣動外形設(shè)計過程中，準(zhǔn)確預(yù)測動穩(wěn)定性導(dǎo)數(shù)非常重要。動穩(wěn)定性導(dǎo)數(shù)可通過工程算法、實驗方法(風(fēng)洞實驗、飛行試驗)以及計算流體力學(xué)(CFD)方法獲得。工程算法只適用于簡單的外形，且精度較低。而實驗方法所需的時間較長，實驗費用比較昂貴，且由干擾因素引起的誤差也是不能避免的。CFD方法作為一種經(jīng)濟且相對準(zhǔn)確的的方法，而被廣泛應(yīng)用于求解飛行器的動穩(wěn)定性導(dǎo)數(shù)。文獻(xiàn)[1-2]通過采用準(zhǔn)定常方法對彈箭穩(wěn)態(tài)錐進運動進行模擬，得到了不同外形彈箭的俯仰組合動導(dǎo)數(shù)。文獻(xiàn)[3-4]采用基于ALE動態(tài)變形網(wǎng)格的非定常方法對帶翼導(dǎo)彈標(biāo)準(zhǔn)模型(BFM)的俯仰和滾轉(zhuǎn)阻尼導(dǎo)數(shù)進行了預(yù)測，仿真結(jié)果與實驗結(jié)果吻合良好。文獻(xiàn)[5-6]使用頻域方法預(yù)測了帶翼導(dǎo)彈在不同條件下的俯仰動導(dǎo)數(shù)，結(jié)果表明當(dāng)縮減頻率在一定范圍內(nèi)取值時，該方法具有和雙時間步方法相當(dāng)?shù)木?，但?dāng)減縮頻率降低到一定值后，結(jié)果不再理想。文獻(xiàn)[7-9]分別采用滑移網(wǎng)格方法和剛性動網(wǎng)格方法對飛機和帶翼導(dǎo)彈的強迫俯仰振動過程進行模擬，得到了模型在不同迎角下的俯仰組合動導(dǎo)數(shù)，該方法可避免網(wǎng)格變形對計算結(jié)果的影響。文獻(xiàn)[10]建立了適用于高速細(xì)長體導(dǎo)彈的滾轉(zhuǎn)動導(dǎo)數(shù)風(fēng)洞實驗裝置，并得到了一致性良好的實驗結(jié)果。文獻(xiàn)[11]結(jié)合CFD仿真方法和Kriging代理模型得到了不同設(shè)計參數(shù)下乘波體的動導(dǎo)數(shù)，并將結(jié)果用于乘波體的荷蘭滾模態(tài)分析。

以上研究對彈箭動導(dǎo)數(shù)求解具有一定借鑒價值，但這些研究只考慮了彈箭的單自由度角運動狀態(tài)，而忽略了彈箭耦合運動帶來的影響。文獻(xiàn)[12-14]的研究結(jié)果表明飛行器耦合運動產(chǎn)生的非定常遲滯特性與單自由度運動產(chǎn)生的遲滯特性不一致，旋轉(zhuǎn)運動對縱向俯仰組合動導(dǎo)數(shù)有顯著影響，使得組合動導(dǎo)數(shù)變得極不穩(wěn)定，同時當(dāng)轉(zhuǎn)速較大時，旋轉(zhuǎn)運動還將顯著增大橫向組合動導(dǎo)數(shù)的非線性。而關(guān)于耦合運動對彈箭動導(dǎo)數(shù)的影響的研究還未見報道。彈箭在實際飛行過程中，通常同時存在俯仰、滾轉(zhuǎn)運動以及彈軸繞速度方向的錐進運動，因此研究彈箭在耦合運動狀態(tài)下的動導(dǎo)數(shù)變化規(guī)律，對研究彈箭的非定常氣動特性和飛行穩(wěn)定性更具有實用價值。

本文在已有的風(fēng)洞實驗數(shù)據(jù)基礎(chǔ)上，提出了一種基于歐拉轉(zhuǎn)動定理和滑移網(wǎng)格技術(shù)的復(fù)雜角運動模擬方法，并將得到的角速度結(jié)果指定給球形滑移網(wǎng)格區(qū)，從而模擬彈箭復(fù)雜角運動過程。通過對得到的氣動力系數(shù)遲滯環(huán)進行參數(shù)辨識，得到了彈箭在耦合運動狀態(tài)下的俯仰組合動導(dǎo)數(shù)和升力系數(shù)動導(dǎo)數(shù)，并分析了滾轉(zhuǎn)頻率和錐進頻率對結(jié)果的影響規(guī)律。

1 動導(dǎo)數(shù)計算方法

基于小幅強迫振動法[7-8]，采用非定常CFD方法模擬滾轉(zhuǎn)運動、錐進運動與強迫俯仰運動耦合的角運動過程，并采用積分法對所得彈箭氣動系數(shù)遲滯環(huán)進行辨識，求解相應(yīng)的彈箭動穩(wěn)定性導(dǎo)數(shù)，以分析耦合運動對其動導(dǎo)數(shù)的影響規(guī)律。

1.1 耦合運動描述

圖1即為本文模擬的彈箭三自由度耦合角運動過程示意圖。圖中坐標(biāo)系OXYZ、OXBYBZB分別為基準(zhǔn)坐標(biāo)系和彈體固連坐標(biāo)系。彈體固連坐標(biāo)系可由基準(zhǔn)坐標(biāo)系經(jīng)過三次旋轉(zhuǎn)得到，即首先繞OX軸逆時針轉(zhuǎn)過φ(進動角)，再繞所得OZ1軸逆時針轉(zhuǎn)α(迎角)，最后繞所得OXB軸逆時針轉(zhuǎn)γ(滾轉(zhuǎn)角)。在如圖所示坐標(biāo)系下，彈箭一邊繞彈軸做勻速滾轉(zhuǎn)運動，一邊在迎角平面內(nèi)做強迫俯仰運動，同時迎角平面繞OX軸勻速偏轉(zhuǎn)。該角運動過程可描述為

其中，f1為彈箭繞OX軸的進動頻率，f2為彈箭做強迫俯仰運動的俯仰頻率，f3為滾轉(zhuǎn)頻率,三個角運動頻率單位均為 Hz。α為彈箭在迎角平面內(nèi)的瞬時迎角，α0為初始迎角，αm為迎角振動幅值，γ為滾轉(zhuǎn)角，φ為進動角。當(dāng)f1、f2和f3取值不同時，彈箭對應(yīng)做不同形式的運動：(1)當(dāng)f2≠0、f1=f3=0時，彈箭做單純的強迫俯仰運動；(2)當(dāng)f1≠0、f2=f3=0時，彈箭做單純的錐進運動；(3)當(dāng)f3≠0、f1=f3=0時，彈箭做單純的滾轉(zhuǎn)運動；(4)當(dāng)f2≠0、f1=0、f3≠0不為0時，彈箭做俯仰滾轉(zhuǎn)耦合角運動；(5)當(dāng)f2≠0、f1≠0、f3=0，彈箭做俯仰錐進耦合運動；(6)當(dāng)f2≠0、f1≠0、f3≠0時，彈箭即做三自由度(俯仰、滾轉(zhuǎn)、錐進)耦合角運動。

圖1 耦合運動示意圖Fig.1 The coupled motion of roll-cone-pitch

由彈體坐標(biāo)系的定義可得，基準(zhǔn)坐標(biāo)系到彈體固連坐標(biāo)系的坐標(biāo)轉(zhuǎn)換矩陣為:

1.2 數(shù)值仿真方法

1.2.1 幾何模型

本文的計算對象為如圖2所示的一種帶有四片尾翼的彈箭。為提高飛行穩(wěn)定性，該彈采用了長尾桿結(jié)構(gòu)，使整個彈的壓心向彈體尾部移動，以提高彈體的穩(wěn)定儲備量。彈箭基本幾何參數(shù)如下：彈徑為D，全彈長為L1=6.8D，彈頭部曲面半徑R=4.57D,尾桿長度為L2=2.6D，尾桿直徑為d=0.48D，尾桿與圓柱部之間通過一個45°的錐形部連接，錐形部長度L3=0.3D;尾翼傾角為0°，翼展為L4=2.8D，翼片寬度為b=0.28D。

圖2 幾何模型Fig.2 Geometry model

1.2.2 控制方程與邊界條件

滑移網(wǎng)格方法是模擬非定常流動過程的常用方法。該方法將流場劃分為外部靜止區(qū)和內(nèi)部滑移區(qū)。計算模型的運動與內(nèi)部滑移區(qū)網(wǎng)格的運動是同步的，因此通過指定內(nèi)部滑移區(qū)網(wǎng)格的運動規(guī)律，即可實現(xiàn)對計算模型運動規(guī)律的模擬。兩區(qū)域在交界面處通過插值計算進行數(shù)據(jù)傳遞?；凭W(wǎng)格方法在計算過程中，不需要進行網(wǎng)格重構(gòu)，計算精度不會受網(wǎng)格變形的影響，具有計算量小易于實現(xiàn)的優(yōu)點。因此，本文采用滑移網(wǎng)格方法對彈箭耦合角運動過程進行模擬。

流場求解采用三維非定常雷諾平均N-S方程，其積分形式如下：

式中V為任意控制體；W是守恒變量；F為無黏通量；Fv為黏性通量；?V為控制體的邊界；n為控制體邊界單位外法向矢量；Re為雷諾數(shù)。

采用有限體積法對方程進行空間離散，其中無黏通量采用Roe進行計算，界面變量采用二階MUSCL格式進行重構(gòu)；黏性通量采用二階中心差分進行計算。湍流模型采用標(biāo)準(zhǔn)k-ε兩方程模型。遠(yuǎn)場來流采用壓力遠(yuǎn)場條件，彈體邊界采用無滑移邊界條件，壁面區(qū)域采用標(biāo)準(zhǔn)壁面函數(shù)，并取彈體直徑D為參考長度，取彈箭橫截面積為參考面積。

求解迭代采用雙時間步推進。每個時間步內(nèi)，迭代次數(shù)取為25次，以保證在每一個時間步內(nèi)計算達(dá)到收斂。每個時間步內(nèi)模型的角速度的計算方法將在1.2.3節(jié)做進一步介紹。前期的計算結(jié)果表明，當(dāng)計算時間達(dá)到3個俯仰運動周期時，俯仰力矩系數(shù)和升力系數(shù)計算結(jié)果均達(dá)到穩(wěn)定振蕩狀態(tài)，因此后續(xù)計算中將第3個俯仰周期的計算結(jié)果作為最終的結(jié)果。

1.2.3 網(wǎng)格劃分

構(gòu)建合理的網(wǎng)格滑移面是滑移網(wǎng)格技術(shù)使用的關(guān)鍵。為滿足模擬彈箭復(fù)雜角運動的需要，將內(nèi)部滑移區(qū)取為球形，外部靜止區(qū)取為圓柱形，兩區(qū)域的交界處的網(wǎng)格滑移面即為球面。將計算模型放置在球心位置，并使計算模型的質(zhì)心與球心重合。這種網(wǎng)格劃分方式的目的在于，當(dāng)球形網(wǎng)格區(qū)域繞球心以任意姿態(tài)轉(zhuǎn)動時，都能保證滑移區(qū)網(wǎng)格和靜止區(qū)網(wǎng)格均不發(fā)生干涉。同時由于模型的質(zhì)心與球心重合，因此內(nèi)部滑移區(qū)網(wǎng)格繞球心的角運動可以直接等同于計算模型繞質(zhì)心的角運動。

采用結(jié)構(gòu)網(wǎng)格對流場進行劃分。遠(yuǎn)場區(qū)域的徑向和軸向高度值分別取為30D和70D，以防止外邊界處反射對結(jié)果的影響。O形網(wǎng)格技術(shù)被用于加密彈體周圍的附面層網(wǎng)格。在黏性邊界層內(nèi)布置了25層網(wǎng)格，第一層網(wǎng)格的厚度被取為10 μm，并且網(wǎng)格延伸率被控制在1.1以下，以保證壁面y+值小于0.5。為進行網(wǎng)格無關(guān)系驗證，按本文提出的網(wǎng)格劃分方案，通過對滑移區(qū)網(wǎng)格進行局部加密得到了網(wǎng)格數(shù)量為200萬、340萬、480萬的三套網(wǎng)格。利用三套網(wǎng)格對來流馬赫數(shù)Ma=0.6，初始迎角為3°，f1=f2=f3=80和f1=80、f2=f3=0兩種非定常運動狀態(tài)進行了仿真求解，結(jié)果表明，當(dāng)網(wǎng)格數(shù)量大于340萬時，所得彈箭各項氣動力系數(shù)的遲滯環(huán)隨網(wǎng)格數(shù)增加無顯著變化，即達(dá)到網(wǎng)格無關(guān)性要求，因此本文采用網(wǎng)格數(shù)為340萬的網(wǎng)格(見圖3)劃分方案進行計算。

圖3 流場網(wǎng)格Fig.3 Grid meshing of flow field

1.2.4 基于歐拉轉(zhuǎn)動定理的角速度計算

采用滑移網(wǎng)格技術(shù)模擬既定的繞心運動，需要指定每個時間步上的角速度矢量ω，若已知時間步長恒定為Δt，則彈箭在第n個時間步的角位移為

式(3)中：η為彈箭第n個時間步的姿態(tài)角矢量，ωi、Δηi分別為第i個時間步內(nèi)的角速度矢量和角位移矢量。

目前，滑移網(wǎng)格技術(shù)在應(yīng)用中角速度矢量ω通常采用求導(dǎo)法計算，通過對既定的運動方程求導(dǎo)并向基準(zhǔn)坐標(biāo)系投影從而得到ω的表達(dá)式。由坐標(biāo)轉(zhuǎn)換關(guān)系可得彈箭角速度矢量ω表達(dá)式：

通過對式(1)求導(dǎo)并結(jié)合坐標(biāo)轉(zhuǎn)換矩陣，將結(jié)果投影到基準(zhǔn)坐標(biāo)系(OXYZ)可得ωd的表達(dá)式為:

上述角速度計算方法相對簡單，但實際仿真過程中發(fā)現(xiàn)，利用求導(dǎo)法求得的角速度進行模擬并不能得到理想的結(jié)果。因為由求導(dǎo)法計算的角速度矢量僅是每個時間步起始點的值，而時間步長Δt為有限值，因此在每個時間步內(nèi)必然存在積分誤差，且時間步長越大誤差越大，且隨著計算周期的增加，彈箭角運動的累積誤差也隨之增加，最終彈箭運動規(guī)律嚴(yán)重偏離式(1)給出的角運動規(guī)律。

為此本文提出了基于歐拉轉(zhuǎn)動定理的角運動計算方法，其實質(zhì)是通過對每個時間步內(nèi)的角速度進行修正，從而消除每個時間步內(nèi)的角位移誤差。推導(dǎo)過程如下：設(shè)第i個時間步開始時刻為t1，由式(1)可得姿態(tài)角為α1、γ1、φ1，將姿態(tài)角帶入式(2)可得轉(zhuǎn)換矩陣為NAB1(α1，γ1，φ1)，并記對應(yīng)的彈體固連坐標(biāo)系為LB1；設(shè)第i+1個時間步開始時刻(即第i個時間步結(jié)束時刻)為t2=t1+Δt，同樣的，由式(1)可得姿態(tài)角為α2、γ2、φ2，將姿態(tài)角帶入式(2)可得對應(yīng)轉(zhuǎn)換矩陣為NAB2(α2，γ2，φ2)，并記錄對應(yīng)的彈體固連坐標(biāo)系為LB2；由此可得從LB1到LB2的坐標(biāo)系轉(zhuǎn)換矩陣表達(dá)式為:

式中，nij為轉(zhuǎn)換矩陣各分量，為節(jié)約篇幅這里不具體給出。

另一方面，由歐拉轉(zhuǎn)動定理可知，彈箭第i+1時間步對應(yīng)的彈體固連坐標(biāo)系，一定可由第i時間步對應(yīng)的彈體固連系經(jīng)過一次旋轉(zhuǎn)得到，即必存在單位矢量u(ux,uy,uz)，以及使LB1繞單位矢量u經(jīng)過一次旋轉(zhuǎn)與LB2重合的角度θ。此時對應(yīng)的轉(zhuǎn)換矩陣為著名的羅德里格斯轉(zhuǎn)換矩陣:

根據(jù)N12=R可聯(lián)立解得:

當(dāng)時間步較小時，近似取第i時間步的角速度矢量為ω=u·(θ/Δt)，并向基準(zhǔn)坐標(biāo)系(OXYZ)投影即可得到由歐拉方法計算的近似角速度矢量:

在利用滑移網(wǎng)格求解過程中，角速度矢量分別采用求導(dǎo)法(式5)和歐拉轉(zhuǎn)動法(式9)在每個時間步為內(nèi)部滑移網(wǎng)格區(qū)指定角速度。對式(1)給出的耦合角運動進行模擬求解，其中取f1=30、f2=40、f3=20，并提取了彈箭在仿真過程中的實際角運動規(guī)律結(jié)果如圖4所示。由圖可知，利用求導(dǎo)法得到彈箭俯仰角與式(1)給出的正弦曲線存在較大偏差，且相對誤差呈周期振蕩；而由歐拉轉(zhuǎn)動法得到的彈箭俯仰角變化規(guī)律，與正弦曲線完全重合。結(jié)果表明，利用該方法可實現(xiàn)對彈箭復(fù)雜角運動的精確模擬，從而提高計算精度。值得注意的是，本文所提出的角運動模擬方法，不僅適用于俯仰滾轉(zhuǎn)耦合運動，還適用于更為復(fù)雜的已知運動方程的繞心運動。

圖4 求導(dǎo)法和歐拉轉(zhuǎn)動法對角運動的模擬精度對比Fig.4 Comparation of simulation accuracy on angle motion between derivative method and Euler rotation method

1.3 動導(dǎo)數(shù)辨識方法

為求解計算彈箭在耦合角運動狀態(tài)下的俯仰組合動導(dǎo)數(shù)和升力系數(shù)動導(dǎo)數(shù)，需要對CFD仿真得到的氣動系數(shù)遲滯環(huán)進行參數(shù)辨識。采用積分法對氣動系數(shù)遲滯環(huán)進行辨識求解，具體推導(dǎo)過程如下:

由式(1)可得彈箭在迎角平面內(nèi)的迎角運動方程:

式中，ωa=2πf2為俯仰角速度圓頻率；q為瞬時俯仰角速度。

彈箭做小幅強迫俯仰振動時，其非定常氣動力和力矩的泰勒展開式為：

將式(10)帶入式(11)，利用三角函數(shù)系的正交性，在一個周期上進行積分，并進行無量綱化處理(詳細(xì)推導(dǎo)過程可參見文獻(xiàn)[15,16]),可得俯仰組合動導(dǎo)數(shù)的計算公式為:

類似的通過對彈箭升力系數(shù)遲滯環(huán)進行積分可得升力系數(shù)組合動導(dǎo)數(shù)C2計算公式為

需要指出的是，仿真直接得到的通常是在基準(zhǔn)坐標(biāo)下的空氣動力和力矩系數(shù)的結(jié)果.本文討論的俯仰力矩系數(shù)和升力系數(shù)均是相對于迎角平面，其中俯仰力矩的方向垂直于迎角平面(彈軸OXB和OX組成的平面)，升力方向位于迎角平面內(nèi)并垂直于OX軸，因此需要將俯仰力矩系數(shù)和升力系數(shù)向迎角平面投影，即:

將式(14)、式(15)分別帶入式(12)、式(13)可得彈箭在迎角平面內(nèi)的俯仰組合動導(dǎo)數(shù)和升力系數(shù)動導(dǎo)數(shù)計算公式分別為:

2 方法驗證

2.1 BFM標(biāo)模算例驗證

采用帶翼導(dǎo)彈標(biāo)準(zhǔn)模型BFM[9]動導(dǎo)數(shù)實驗數(shù)據(jù)檢驗本文所采用的CFD仿真方法以及動導(dǎo)數(shù)辨識方法的準(zhǔn)確性。基于本文1.2節(jié)所述數(shù)值仿真方法以及1.3節(jié)所述動導(dǎo)數(shù)辨識方法，對BFM標(biāo)模在馬赫數(shù)為1.58、1.76、1.89、2.16、2.48，起始迎角為1.5°、振幅為1.5°條件下的強迫俯仰振動進行了求解計算(仿真過程中取f1=f3=0，f2≠0，即可模擬彈箭的強迫俯仰振動)。參考長度取BFM標(biāo)模的圓柱部直徑D0，參考面積取圓柱部橫截面積，質(zhì)心位置到彈頭部距離為5D0。將俯仰組合動導(dǎo)數(shù)的計算結(jié)果與文獻(xiàn)[17,18]的實驗結(jié)果進行對比，結(jié)果如圖5所示。由圖5可知，俯仰組合動導(dǎo)數(shù)計算結(jié)果與實驗結(jié)果變化規(guī)律一致，兩者絕對值均隨馬赫數(shù)增加而減小，其中仿真值略小于實驗值，相對誤差在Ma=1.89時最大為6.4%。結(jié)果表明，所采用數(shù)值仿真方法和動導(dǎo)數(shù)辨識方法對帶尾翼彈箭的動導(dǎo)數(shù)求解具有較高精度。

圖5 BFM標(biāo)模俯仰組合動導(dǎo)數(shù)隨馬赫數(shù)變化曲線Fig.5 Curve of pitching and rolling damping derivative of BFM model with Mach number

2.2 數(shù)值方法實驗驗證

由于缺少計算模型的動導(dǎo)數(shù)實驗數(shù)據(jù)，無法直接對該模型的動導(dǎo)數(shù)計算精度進行分析。作為補充，采用已有的滾轉(zhuǎn)風(fēng)洞實驗數(shù)據(jù)對非定常數(shù)值方法的計算精度進行檢驗。實驗?zāi)Ｐ蛶缀谓Y(jié)構(gòu)與圖2一致，但尾翼傾角為13°(圖2中的模型尾翼傾角為0°)。風(fēng)洞實驗在南京理工大學(xué)HG-4號風(fēng)洞進行，實驗安裝圖如圖6所示。彈體通過隔轉(zhuǎn)軸承安裝在天平軸上，使彈體可繞天平軸進行自由滾轉(zhuǎn)。實驗過程中，實驗?zāi)Ｐ驮谛敝梦惨懋a(chǎn)生的導(dǎo)旋力矩下做滾轉(zhuǎn)運動，并最終達(dá)到穩(wěn)定轉(zhuǎn)速。實驗測量了Ma=0.6時模型在不同迎角下的穩(wěn)定轉(zhuǎn)速和所受空氣動力，表1給出了Ma=0.6時模型在不同迎角下所達(dá)到的穩(wěn)定轉(zhuǎn)速。

圖6 風(fēng)洞實驗安裝圖Fig.6 Wind tunnel test setup

采用1.2節(jié)所述的非定?；凭W(wǎng)格方法模擬彈箭的滾轉(zhuǎn)運動過程，在仿真過程中，取f1=f2=0，f3按表1給出的結(jié)果取值，以保證仿真條件與實驗條件的一致。升力系數(shù)仿真結(jié)果和實驗結(jié)果進行對比，結(jié)果如圖7所示。由圖可知，數(shù)值仿真結(jié)果與實驗結(jié)果吻合較好，最大相對誤差小于10%，表明本文采用的數(shù)值仿真方法能對計算模型的氣動參數(shù)進行準(zhǔn)確預(yù)測。

圖7 Ma=0.6時，升力系數(shù)隨迎角變化曲線Fig.7 Curve of lift coefficient with the angle of attack (Ma=0.6)

迎角/(°)-2048平衡轉(zhuǎn)速/Hz216212210190

3 結(jié)果與分析

采用本文1.2節(jié)所述角運動模擬方法，對彈箭滾轉(zhuǎn)、錐進運動與強迫俯仰運動相耦合的角運動過程進行模擬。計算條件如下：取計算馬赫數(shù)Ma為0.65，起始迎角α0取為3°，迎角振蕩幅值αm取為1.5°，質(zhì)心位置到彈頂距離取為0.4D。仿真過程中，彈箭俯仰振蕩頻率f2取為定值20 Hz，錐進頻率f1和滾轉(zhuǎn)頻率f3取值范圍為0到100 Hz，并制定了如表2所示的4組仿真方案，記為F1、F2、F3、F4，其中F1方案用于分析彈箭在俯仰滾轉(zhuǎn)耦合運動狀態(tài)下滾轉(zhuǎn)頻率對動導(dǎo)數(shù)辨識結(jié)果的影響，F(xiàn)3方案用于分析彈箭在錐進俯仰耦合狀態(tài)下錐進頻率對動導(dǎo)數(shù)辨識結(jié)果的影響，F(xiàn)2和F4方案用于分析彈箭同時存在錐進、滾轉(zhuǎn)以及俯仰運動的三自由度耦合運動狀態(tài)下錐進頻率和滾轉(zhuǎn)頻率對動導(dǎo)數(shù)辨識結(jié)果的影響。

表2 仿真方案Table 2 Simulation scheme

3.1 俯仰力矩系數(shù)和升力系數(shù)遲滯環(huán)

從四組仿真方案中選取了三種具有代表性的情況的結(jié)果，用于分析滾轉(zhuǎn)運動和錐進運動對彈箭氣動力遲滯環(huán)以及繞流特性的影響，分別記為P1、P2、P3、P4：① P1是僅存在強迫俯仰振蕩的情況(f1=0 Hz，f3=0 Hz，f2=20 Hz)；② P2是滾轉(zhuǎn)運動與強迫俯仰振蕩耦合的情況(f1=0 Hz，f3=80 Hz，f2=20 Hz)；③ P3是錐進運動與強迫俯仰振蕩耦合的情況(f1=80 Hz，f3=0 Hz，f2=20 Hz)。④ 補充了同時存在俯仰滾轉(zhuǎn)和錐進運動(f1=80 Hz，f3=80 Hz，f2=20 Hz)的情況進行對比分析，并記為P4。

圖8給出了四種條件下彈箭俯仰力矩系數(shù)滯回曲線仿真結(jié)果。由圖可知，P2、P3、P4的結(jié)果的絕對值整體小于P1的結(jié)果，其中P4尤為顯著，說明在彈箭做強迫俯仰振動過程中，滾轉(zhuǎn)運動或錐進運動的存在，均會導(dǎo)致彈箭俯仰力矩系數(shù)絕對值減小，當(dāng)同時存在滾轉(zhuǎn)運動和錐進運動時，變化尤為顯著。P2、P3、P4的結(jié)果與P1的結(jié)果相比，當(dāng)存在滾轉(zhuǎn)運動或錐進運動時，俯仰力矩系數(shù)滯回曲線不再光滑，而是出現(xiàn)了明顯的振蕩和偏移。其中由滾轉(zhuǎn)運動(P2)引起的遲滯環(huán)振蕩相比于由錐進運動(P3)引起的振蕩更為顯著。這是因為在彈箭在強迫俯仰運動過程中，滾轉(zhuǎn)運動或錐進運動的存在均會引起彈身及尾翼周圍的流場不再對稱，產(chǎn)生周期性的干擾力和力矩，使遲滯環(huán)發(fā)生振蕩。滾轉(zhuǎn)運動還會引起尾翼相對于迎角平面的方位角不斷改變，這將進一步增強俯仰力矩系數(shù)的振蕩。而本文所定義的錐進運動，并不會引起尾翼相對于迎角平面的方位角發(fā)生改變，因此錐進運動引起的遲滯環(huán)振蕩相對較弱。從圖中還可以看出，P3的結(jié)果的絕對值小于P2，表明在滾轉(zhuǎn)頻率和錐進頻率取值相同的情況下，由錐進運動引起的俯仰力矩系數(shù)的減小更為顯著。此外，P4的遲滯環(huán)的振蕩和偏移都非常顯著，說明錐進運動和滾轉(zhuǎn)運動對結(jié)果的影響存在疊加效果。圖9為四種條件下彈箭升力系數(shù)滯回曲線仿真結(jié)果。將圖9和圖8對比可知，升力系數(shù)滯回曲線的變化規(guī)律和俯仰力矩系數(shù)滯回曲線的變化規(guī)律是相似的，只是在數(shù)值上有一定差異。滾轉(zhuǎn)運動或錐進運動的存在均會導(dǎo)致彈箭升力數(shù)遲滯出現(xiàn)明顯的振蕩和偏移，且當(dāng)滾轉(zhuǎn)和錐進運動同時存在時，這種變化尤為顯著。以上分析表明，彈箭在耦合角運動狀態(tài)下的氣動力和力矩系數(shù)遲滯環(huán)，與單純的強迫俯仰振動得到的遲滯環(huán)相比存在顯著差異，這必然對彈箭動導(dǎo)數(shù)的穩(wěn)定性辨識結(jié)果產(chǎn)生影響。

圖8 俯仰力矩系數(shù)滯回曲線Fig.8 Hysteresis curve of pitching moment coefficient

3.2 動導(dǎo)數(shù)辨識結(jié)果

為定量分析彈箭滾轉(zhuǎn)運動對其俯仰動導(dǎo)數(shù)辨識結(jié)果的影響程度，采用1.3節(jié)所述方法對俯仰力矩系數(shù)和升力系數(shù)遲滯曲線進行參數(shù)辨識，得到結(jié)果如圖10～圖13所示。

圖10和圖11分別是根據(jù)表1給出的仿真方案得到的俯仰組合動導(dǎo)數(shù)隨滾轉(zhuǎn)頻率和錐進頻率的變化曲線。由圖可知，俯仰組合動導(dǎo)數(shù)隨滾轉(zhuǎn)頻率和錐進頻率增加均具有減小趨勢，且滾轉(zhuǎn)頻率引起的減小幅值更為顯著。F2和F4的結(jié)果均分別小于F1和F3的結(jié)果，表明由滾轉(zhuǎn)運動引起的結(jié)果的減小與由錐進運動引起的結(jié)果的減小，存在疊加效果。但從圖中可以看出，F(xiàn)2與F1之間的差值以及F4與F3之間的差值分別隨滾轉(zhuǎn)頻率增加而有所增大，其中圖11尤為明顯，表明當(dāng)同時存在滾轉(zhuǎn)運動和錐進運動時，由滾轉(zhuǎn)運動引起的結(jié)果的減小與由錐進運動引起的結(jié)果的減小，并非簡單的線性疊加關(guān)系，即滾轉(zhuǎn)運動和錐進運動對彈箭俯仰動導(dǎo)數(shù)的影響存在耦合效應(yīng)。

圖10 俯仰組合動導(dǎo)數(shù)隨滾轉(zhuǎn)頻率的變化曲線Fig.10 Pitch damping derivative changing with the roll frequency

圖11 俯仰組合動導(dǎo)數(shù)隨錐進頻率的變化曲線Fig.11 Pitch damping derivative changing with the coning frequency

圖12和圖13分別是根據(jù)表1給出的仿真方案得到的升力系數(shù)動導(dǎo)數(shù)隨滾轉(zhuǎn)頻率和錐進頻率的變化曲線。由圖可知，F(xiàn)2和F4的結(jié)果均分別小于F1和F3的結(jié)果，表明滾轉(zhuǎn)運動和錐進運動對升力系數(shù)的影響同樣存在的疊加效果。圖12的結(jié)果顯示，彈箭升力系數(shù)動導(dǎo)數(shù)隨滾轉(zhuǎn)頻率增加而明顯的減小趨勢。其中僅存在滾轉(zhuǎn)運動時(F1)，升力系數(shù)動導(dǎo)數(shù)隨滾轉(zhuǎn)頻率增加而單調(diào)減小，且變化規(guī)律是接近線性的。當(dāng)同時存在錐進運動時(F2)，升力系數(shù)變化規(guī)律出現(xiàn)振蕩，這是因為滾轉(zhuǎn)運動和錐進運動與強迫俯仰運動的共同作用使彈箭繞流特性更加復(fù)雜，尤其在尾錐以及尾翼等部位的壓力分布出現(xiàn)了非線性變化，從而引起了彈箭整體所受空氣動力的改變。由圖13可知，升力系數(shù)動導(dǎo)數(shù)隨錐進頻率增加不再單調(diào)減小，而是按先減小后增大的規(guī)律變化，變化規(guī)律具有明顯非線性特性，且當(dāng)錐進頻率F1=40 Hz時，升力系數(shù)動導(dǎo)數(shù)達(dá)到最小值。圖12和圖13結(jié)果的差異，體現(xiàn)了滾轉(zhuǎn)運動與錐進運動對動導(dǎo)數(shù)辨識結(jié)果的影響存在明顯的差異。

圖12 升力系數(shù)動導(dǎo)數(shù)隨滾轉(zhuǎn)頻率的變化曲線Fig.12 Lift coefficient derivative changing with the roll frequency

圖13 俯仰組合動導(dǎo)數(shù)隨錐進頻率的變化曲線Fig.13 Lift coefficient derivative changing with the coning frequency

表3給出了滾轉(zhuǎn)力矩系數(shù)和升力系數(shù)的部分計算結(jié)果，表中的相對變化率是相對于僅存在強迫俯仰運動(f1=0 Hz，f2=0 Hz)的情況而言。從表中可以看出，除f1=80 Hz，f2=0 Hz(即錐進運動與強迫俯仰運動耦合)的情況外，俯仰組合動導(dǎo)數(shù)和升力系數(shù)動導(dǎo)數(shù)均達(dá)到10%以上，其中當(dāng)滾轉(zhuǎn)頻率f3=80 Hz、錐進頻率f1=40 Hz時，俯仰組合動導(dǎo)數(shù)減小了23.1%，升力系數(shù)動導(dǎo)數(shù)減小了21.2%，結(jié)果表明，彈箭在三自由度耦合角運動(滾轉(zhuǎn)運動和錐進運動與強迫俯仰運動相耦合)狀態(tài)下，動穩(wěn)定性導(dǎo)數(shù)的辨識結(jié)果與單純的強迫俯仰振動的結(jié)果相比具有顯著差異，因此在進行彈箭動導(dǎo)數(shù)計算和穩(wěn)定性分析時需被充分考慮。

表3 不同條件下的動導(dǎo)數(shù)計算結(jié)果對比Table 3 Comparation of dynamic derivative under different conditions

3.3 壓力分布

為進一步分析滾轉(zhuǎn)運動和錐進運動對彈箭所受空氣動力的影響，提取了四種情況(P1～P4)下，彈箭尾錐位置和水平尾翼位置的壓力分布情況，結(jié)果如圖14和圖15所示。圖14是彈箭在瞬時迎角為3°尾錐部的周向壓力分布情況。由圖14可知，當(dāng)彈箭只作單純的強迫俯仰振動時(無滾轉(zhuǎn)無錐進，對應(yīng)曲線P1)，彈箭尾錐部左右兩側(cè)壓力對稱，且迎風(fēng)面壓力明顯大于背風(fēng)面。當(dāng)彈箭存在滾轉(zhuǎn)運動時(P2)或錐進運動時(P3)，尾錐部壓力分布不再對稱，主要表現(xiàn)為彈箭右側(cè)壓力大于左側(cè)壓力，且迎風(fēng)面和背風(fēng)面壓力差亦發(fā)生改變。原因為彈箭尾錐部錐度較大，彈箭繞流在尾錐部形成明顯的渦。當(dāng)彈箭無滾轉(zhuǎn)無錐進時，尾錐部的渦較對稱，此時尾錐部的壓力分布是對稱的；當(dāng)彈箭存在滾轉(zhuǎn)或錐進運動時，由于空氣黏性的作用，彈箭附面層位移厚度和渦結(jié)構(gòu)將發(fā)生非對稱畸變，即產(chǎn)生馬格努斯效應(yīng)，使尾錐部周向壓力分布發(fā)生改變。此外，對比P2和P3的結(jié)果可知，錐進運動產(chǎn)生的非對稱性更為顯著。當(dāng)彈箭同時存在滾轉(zhuǎn)和錐進運動時(P4)，尾錐部壓力分布非對稱性進一步增加，表明滾轉(zhuǎn)運動和錐進運動對彈箭繞流特性的影響存在耦合效應(yīng)，且這種耦合效應(yīng)是非線性的。

圖14 尾錐部的壓力分布Fig.14 Pressure distribution of the tail cone

圖15給出了彈箭在瞬時迎角為3°時，水平尾翼的壓力分布情況(考慮到升力系數(shù)和俯仰力矩系數(shù)主要受水平尾翼的壓力分布的影響，因此只提取了水平尾翼的壓力分布)。由圖可知，當(dāng)無滾轉(zhuǎn)且無錐進時(P1)，由于迎角存在，左右尾翼的下側(cè)翼面壓力均大于上側(cè)壓力，且壓力分布是左右對稱的。當(dāng)存在滾轉(zhuǎn)運動(P2)或錐進運動(P3)時，尾翼翼面微元的瞬時迎角發(fā)生改變，并與初始迎角疊加，使尾翼壓力分布不再對稱，表現(xiàn)為右側(cè)翼面壓力值增大，左側(cè)壓力系數(shù)由正值變?yōu)樨?fù)值。對比曲線P2和曲線P3可知，左側(cè)翼在滾轉(zhuǎn)頻率f3=80 Hz和錐進頻率f1=80 Hz兩種情況下壓力系數(shù)是接近相等的，而右側(cè)翼在兩種情況下的壓力系數(shù)顯著不同。表明由滾轉(zhuǎn)運動和錐進運動引起的水平尾翼壓力分布變化是不同的。當(dāng)同時存在滾轉(zhuǎn)和錐進時(P4)，左右兩翼片的壓力系數(shù)絕對值均大于僅存在滾轉(zhuǎn)或錐進運動的情況。表明滾轉(zhuǎn)運動和錐進運動對尾翼壓力分布的影響存在耦合效應(yīng)，且這種耦合效應(yīng)是非線性的。尾翼壓力分布的改變必然引起彈箭升力系數(shù)以及俯仰力矩系數(shù)發(fā)生改變。

圖15 水平尾翼壓力分布Fig.15 Pressure distribution of the horizontal tail

以上分析表明，當(dāng)存在一定迎角時，滾轉(zhuǎn)運動和錐進運動，使彈箭表面(尤其在彈箭尾錐和尾翼處)的壓力分布不再對稱，即產(chǎn)生馬格努斯效應(yīng)，這必然引起彈箭所受空氣動力發(fā)生改變。由于彈箭在耦合角運動狀態(tài)下，彈箭迎角隨強迫俯仰運動過程不斷改變，由滾轉(zhuǎn)運動和錐進運動引起的馬格努斯效應(yīng)也將隨之改變，從而形成周期性的誘導(dǎo)力和力矩。此外，滾轉(zhuǎn)運動還會引起尾翼相對于迎角平面的相位角不斷改變，引起彈箭所受空氣動力發(fā)生振蕩。以上兩方面的原因?qū)е聫椉隈詈辖沁\動狀態(tài)下的氣動力系數(shù)遲滯環(huán)與單純的強迫俯仰過程相比發(fā)生了明顯振蕩和偏移，從而對彈箭動穩(wěn)定性導(dǎo)數(shù)的辨識結(jié)果產(chǎn)生影響。

4 結(jié) 論

通過對彈箭在俯仰滾轉(zhuǎn)耦合角運動狀態(tài)下的氣動特性進行數(shù)值仿真研究，得出以下結(jié)論：

1) 本文提出的基于歐拉轉(zhuǎn)動定理和滑移網(wǎng)格技術(shù)的復(fù)雜角運動模擬方法，可在彈箭非定常流場求解過程中，有效消除彈箭姿態(tài)角計算的累積誤差。此方法不僅適用于模擬彈箭俯仰滾轉(zhuǎn)耦合角運動，而且對任意給定形式的角運動同樣適用。將數(shù)值仿真結(jié)果與已有的實驗數(shù)據(jù)進行對比驗證了所采用的數(shù)值仿真方法以及動導(dǎo)數(shù)辨識方法具有較高的準(zhǔn)確性；

2) 彈箭在三自由度耦合角運動(滾轉(zhuǎn)運動和錐進運動與強迫俯仰運動相耦合)狀態(tài)下，氣動力系數(shù)遲滯環(huán)發(fā)生明顯的振蕩和偏移，且動穩(wěn)定性導(dǎo)數(shù)的辨識結(jié)果與單純的強迫俯仰振動的結(jié)果相比具有顯著差異，其中，當(dāng)滾轉(zhuǎn)頻率f3=80、錐進頻率f1=40時，俯仰組合動導(dǎo)數(shù)減小了23.1%，升力系數(shù)動導(dǎo)數(shù)減小了21.2%，因此在進行彈箭動導(dǎo)數(shù)計算和穩(wěn)定性分析時需被充分考慮。