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

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

0 引 言

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

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

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

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

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

1.1 耦合運動描述

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

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

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

由彈體坐標系的定義可得，基準坐標系到彈體固連坐標系的坐標轉換矩陣為:

1.2 數(shù)值仿真方法

1.2.1 幾何模型

本文的計算對象為如圖2所示的一種帶有四片尾翼的彈箭。為提高飛行穩(wěn)定性，該彈采用了長尾桿結構，使整個彈的壓心向彈體尾部移動，以提高彈體的穩(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ū)和內部滑移區(qū)。計算模型的運動與內部滑移區(qū)網(wǎng)格的運動是同步的，因此通過指定內部滑移區(qū)網(wǎng)格的運動規(guī)律，即可實現(xiàn)對計算模型運動規(guī)律的模擬。兩區(qū)域在交界面處通過插值計算進行數(shù)據(jù)傳遞?；凭W(wǎng)格方法在計算過程中，不需要進行網(wǎng)格重構，計算精度不會受網(wǎng)格變形的影響，具有計算量小易于實現(xiàn)的優(yōu)點。因此，本文采用滑移網(wǎng)格方法對彈箭耦合角運動過程進行模擬。

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

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

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

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

1.2.3 網(wǎng)格劃分

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

采用結構網(wǎng)格對流場進行劃分。遠場區(qū)域的徑向和軸向高度值分別取為30D和70D，以防止外邊界處反射對結果的影響。O形網(wǎng)格技術被用于加密彈體周圍的附面層網(wǎng)格。在黏性邊界層內布置了25層網(wǎng)格，第一層網(wǎng)格的厚度被取為10 μm，并且網(wǎng)格延伸率被控制在1.1以下，以保證壁面y+值小于0.5。為進行網(wǎng)格無關系驗證，按本文提出的網(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)進行了仿真求解，結果表明，當網(wǎng)格數(shù)量大于340萬時，所得彈箭各項氣動力系數(shù)的遲滯環(huán)隨網(wǎng)格數(shù)增加無顯著變化，即達到網(wǎng)格無關性要求，因此本文采用網(wǎng)格數(shù)為340萬的網(wǎng)格(見圖3)劃分方案進行計算。

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

1.2.4 基于歐拉轉動定理的角速度計算

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

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

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

通過對式(1)求導并結合坐標轉換矩陣，將結果投影到基準坐標系(OXYZ)可得ωd的表達式為:

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

2 方法驗證

2.1 BFM標模算例驗證

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

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

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

由于缺少計算模型的動導數(shù)實驗數(shù)據(jù)，無法直接對該模型的動導數(shù)計算精度進行分析。作為補充，采用已有的滾轉風洞實驗數(shù)據(jù)對非定常數(shù)值方法的計算精度進行檢驗。實驗模型幾何結構與圖2一致，但尾翼傾角為13°(圖2中的模型尾翼傾角為0°)。風洞實驗在南京理工大學HG-4號風洞進行，實驗安裝圖如圖6所示。彈體通過隔轉軸承安裝在天平軸上，使彈體可繞天平軸進行自由滾轉。實驗過程中，實驗模型在斜置尾翼產生的導旋力矩下做滾轉運動，并最終達到穩(wěn)定轉速。實驗測量了Ma=0.6時模型在不同迎角下的穩(wěn)定轉速和所受空氣動力，表1給出了Ma=0.6時模型在不同迎角下所達到的穩(wěn)定轉速。

圖6 風洞實驗安裝圖Fig.6 Wind tunnel test setup

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

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

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

3 結果與分析

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

表2 仿真方案Table 2 Simulation scheme

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

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

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

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

3.2 動導數(shù)辨識結果

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

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

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

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

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

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

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

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

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

3.3 壓力分布

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

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

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

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

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

4 結 論

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

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

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