俯仰阻尼導數分量的CFD數值模擬

席柯，閻超*，黃宇，王文，袁武

(1.北京航空航天大學 航空科學與工程學院，北京100191;2.中國科學院 計算機網絡信息中心，北京100190)

飛行器動態穩定性導數(工程上常稱為“動導數”)是飛行器動態氣動特性設計、彈道設計及動態穩定性分析中的關鍵參數.動導數是分析飛行器動態品質和設計操控系統的原始氣動數據，動導數的準確評估對飛行器設計和飛行都有重要的意義.

傳統上，采用繞定軸俯仰振動(強迫振動或自由振動)方法得到的是俯仰阻尼導數的組合形式Cmq+，即直接阻尼導數(或稱旋轉導數)Cmq和洗流時差導數(或稱加速度導數)的組合.傳統位勢流理論認為:時差導數在組合導數中所占的比例較小，采用Cmq+代替Cmq誤差不大.但是，特別是隨著現代科技的進步和戰爭需求的升級，出現了具有大迎角、超機動飛行能力的飛行器，其外形設計及運動方式都比傳統飛行器復雜，流動的非線性及非定常性也更強.上述結論不一定成立，比如大升力體飛行器的時差導數在組合導數中所占比例較大［1］(可達40% ～50%)，并且實際飛行中，飛行器的俯仰角速度和迎角變化率并不總是相等，因此需要分別確定組合導數的各個分量.

目前風洞試驗及計算對俯仰阻尼導數分量的研究并不多.一般風洞提供的是直勻流，不能進行直接阻尼導數Cmq的試驗.國內，任玉新［2］求解流動控制方程的敏感性方程，直接得到阻尼導數的各個分量，并研究了減縮頻率的影響，劉偉等人［3］采用強迫沉浮運動形式，數值模擬了高超聲速返回艙及彈道外形(HBS)的加速度導數，研究指出在組合導數小于零的情況下，加速度導數的符號可能大于零，起負阻尼作用，屬于動不穩定情況.孫海生［4］在低速風洞中進行了加速度導數的試驗，指出在α＜30°范圍內，加速度導數大約是組合導數的30% ～50%.國外，美國ARL(Army Research Laboratory)的 Weinacht［5-6］及英國的 Qin等人［7］采用定常方法求解錐轉運動及螺旋運動，獲得了ANSR(Army-Navy Spinner Rocket)和裙狀發射體的直接阻尼導數Cmq及時差導數，但該方法只適用于軸對稱外形.

本文結合計算流體力學(CFD)技術開展了飛行器俯仰阻尼導數各個分量的數值模擬研究.采用準定常方法直接獲得直接阻尼導數Cmq，求解強迫沉浮運動獲得洗流時差導數，求解強迫俯仰振動獲得俯仰阻尼組合導數Cmq+.通過HBS及基本帶翼導彈(Basic Finner)標模進行了數值方法驗證及研究，探討了質心位置對各阻尼導數的影響以及大攻角下各阻尼導數的變化規律，在此基礎上計算了日本的再入飛行器Hyflex的各阻尼導數.

1 計算方法

1.1 流動控制方程

流動控制方程為三維非定常可壓縮Navier-Stokes方程.一般曲線坐標系中，無量綱化的方程守恒形式為［8］

式中，Q為守恒變矢量;F，G和H為對流通量;Fv，Gv，Hv為黏性通量;ξ，η 和 ζ為 3 個貼體坐標系方向;t為時間;Re∞為自由來流雷諾數.

流場求解采用基于結構網格的有限體積法，空間格式采用Roe格式和Minmod限制器，時間離散方法為LU-SGS(Lower-Upper Symmetric Gauss-Seidel)方法，非定常時間推進采用雙時間步方法.采用SA湍流模型模擬大攻角下可能存在的分離流動，使用Open-MP并行技術提高計算效率.

1.2 直接阻尼導數計算方法

直接阻尼導數的物理意義是由于飛行器的姿態角速度的變化引起各處局部迎角的變化，結果引起質心前后升力之差而產生附加的阻尼力矩，該阻尼力矩對角速度的導數即為直接阻尼導數，它本質上是準定常量.

設體軸坐標系x軸朝前，y軸朝上，z軸滿足右手法則.本文采用準定常方法，假設飛行器繞z軸以恒定俯仰角速度q拉起，某瞬時攻角α下有

式中，Cm為飛行器以勻角速率q拉起時瞬時攻角α處的俯仰力矩;Cm0為俯仰角速度為零時攻角α處的俯仰力矩;l為特征長度;u∞為來流速度;

Cm0通過定常計算獲得;Cm的計算需要考慮俯仰角速度引起的牽連速度的影響以反映局部迎角的變化，本文在空間格式及壁面邊界條件中加入牽連速度的貢獻，同時網格保持不動，采用準定常方法進行時間推進，經過一段暫態效應，最終力和力矩收斂到定常結果.

需要說明的是，本文采用的方法不同于文獻［5-7］中采用的方法，文獻中Weinacht等人在非慣性坐標系下通過在流動控制方程中引入源項來考慮牽連速度的影響，而本文方法是在慣性坐標系下，求解帶運動速度物體的瞬態情況，只需將現有求解強迫運動的程序稍作修改即可實現.

1.3 洗流時差導數計算方法

洗流時差導數本質上是非定常量，反映了洗流時差對飛行器的阻尼特性.洗流時差導數的作用體現在飛行器受陣風或直接力作用時飛行器動態穩定性的變化，它是飛行器設計時一個重要的設計參數.

在體軸坐標系，給定沉浮運動形式如下［3］:

在式(4)描述的運動形式下，可得瞬時攻角的運動形式為

式中，α0為起始攻角;αm為攻角振幅;θ為俯仰角;k為減縮頻率.

由式(4)、式(5)確定的沉浮運動中，確定運動的狀態變量只有攻角及其各階導數，根據Etkin非定常氣動力模型，Cm可寫成

在小振幅運動情況下，將式(6)泰勒展開并略去高階小量得

俯仰阻尼導數Cmq+的計算采用強迫振動方法，具體數值計算方法參見文獻［9-10］.

2 數值模擬結果及分析

本節采用上述數值計算方法，對HBS標模外形及基本帶翼導彈外形進行計算驗證和研究分析，并在此基礎上對Hyflex飛行器進行研究.

2.1 HBS 外形研究

HBS為一個半球柱、帶有兩段擴張裙部的高超聲速彈道外形標模，如圖1所示.其動態特性有較為精確的試驗結果［11］，常被用來驗證計算結果的準確程度.

圖1 彈道外形示意圖Fig.1 Schematic of hyperballistic shape(HBS)

計算條件為:Ma=6.85，以頭部直徑為參考長度的 Red=0.72 ×106，質心位置 Xcg/L=0.72，α0=0°，5°，10°，15°，20°.強迫沉浮運動及角振動的減縮頻率k定義為:k=ωd/2u∞，d為頭部直徑.本節計算時取 k=0.01，振幅取 αm=1°.

圖2給出了計算得到的俯仰力矩系數曲線.其中圖2(a)為HBS外形俯仰角速度取q=0，100，300，500，700，1 000(°)/s 時的俯仰力矩系數變化曲線，可以看出，俯仰力矩系數隨俯仰角速度呈完全線性變化，其斜率即為直接阻尼導數，斜率都為負，說明直接阻尼導數為負，起正阻尼作用.圖2(b)為強迫沉浮運動時俯仰力矩系數的遲滯環，圖2(c)為強迫角振動時俯仰力矩系數的遲滯環，對比兩圖可以發現，沉浮運動的遲滯環比較瘦，并且0°～15°的遲滯環為順時針方向旋轉，由于遲滯環的旋轉方向決定動導數的正負，面積代表動導數的量值大小，這說明HBS的洗流時差導數在俯仰阻尼導數中所占比例較小并且0°～15°的洗流時差導數為正值，其洗流特性是負阻尼的，屬于動不穩定情況.而相應的強迫振動的遲滯環比較飽滿且各攻角下均為逆時針旋轉.

圖3為計算所得俯仰阻尼導數與試驗值［11］對比，表1為各攻角下直接阻尼導數、洗流時差導數及俯仰阻尼導數的計算結果，其中Cmq+Cmα·為Cmq和直接相加的結果為通過強迫振動直接得到的俯仰阻尼導數，兩者誤差不超過2%，說明本文發展的計算方法是正確的，并且計算結果是精確的.

圖2 彈道外形計算結果Fig.2 Calculation results of hyperballistic shape(HBS)

圖3 彈道俯仰阻尼導數計算結果與試驗結果對比Fig.3 Calculation and experiment results of pitch-damping derivatives for hyperballistic shape(HBS)

表1 彈道阻尼導數計算結果Table1 Caculation results of pitch-damping derivatives for hyperballistic shape(HBS)

2.2 基本帶翼導彈外形研究

基本帶翼導彈為一個尖錐形頭部、圓柱形彈身并帶有4個矩形尾翼的外形，如圖4所示.

圖4 基本帶翼導彈外形Fig.4 Schematic of Basic Finner

計算條件為:Ma=1.96，以底部直徑為參考長度的ReD=0.187×106.強迫沉浮運動及角振動的減縮頻率k定義為:k=ωD/2u∞，D為底部直徑.本節計算時取 k=0.01，振幅取 αm=1°.

圖5展示了固定質心位置Xcg/D=6.1時，導彈的直接阻尼導數Cmq、洗流時差導數及俯仰阻尼導數Cmq+隨攻角的變化規律(圖中曲線采用樣條曲線擬合而成，下同).計算攻角為α=0°，2.5°，4°，5°，10°，15°，20°，25°，30°.本文計算的雷諾數與風洞試驗［12］中的高雷諾數狀態相同.從圖5中可以看出，本文計算的俯仰阻尼導數在大攻角時(10°以上)與試驗值吻合較好，而小攻角時有一定的偏離.因為試驗采用尾支撐方式，Uselton等人［13］注意到小攻角時(6°以下)有較強的支架干擾而導致試驗結果重復度較差，直到攻角增大到一定程度才會減弱支架干擾的影響.直接阻尼導數一直為負值，且其絕對值隨攻角增大而近似線性增大;時差導數在0°攻角附近為正值，在4°攻角處過零點，然后向負向增大，在10°攻角附近達到負向最大，然后緩慢回落，并一直保持負值.時差導數的非線性變化導致了俯仰阻尼導數的非線性變化，但總體來看，導彈的縱向動穩定性隨攻角增大而增強，時差導數符號會有變化，但量值在組合導數中所占比例較小(10%以下).

圖5 各阻尼導數隨攻角變化曲線Fig.5 Variations of pitch-damping derivatives with angles of attack

圖6 各阻尼導數隨質心位置變化曲線Fig.6 Variations of pitch-damping derivatives with center of gravity(CG)location

圖6為0°攻角時計算的導彈各阻尼導數隨質心位置的變化曲線圖.計算的質心位置為Xcg/D=5，6.1，7.從圖6中可以看出，時差導數為正值，在質心Xcg/D=5處幾乎為零且隨質心位置后移而線性增加;直接阻尼導數及組合導數則呈拋物分布，這與文獻［14］給出的質心平移關系式(8)～式(10)吻合很好.

式中scg=(－xcg)/D，坐標定義在體軸系.上式只適用于軸對稱物體.

雖然沒有在此列出，但計算中CNα，CNq，均為正值，因此時差導數數值隨質心后移增加，即朝動不穩定性增大的方向發展.而直接阻尼導數與質心偏移量是復雜的二次函數關系.從數值上來看，存在一個質心位置使俯仰阻尼導數值最大，即動穩定性最差.

2.3 Hyflex外形研究

Hyflex(Hypersonic Flight Experiment)是日本HOPE-X計劃中有關大氣層再入項目的帶翼升力體外形高超聲速飛行器，用于驗證制導和控制以及熱防護材料和結構等技術.它于1996年2月進行飛行試驗，本文選取其飛行末端彈道點(馬赫2.0，高度30 km)進行研究［15］，該彈道點處于飛行器減速傘開傘前，其動態特性對于減速傘的開啟有著至關重要的影響.飛行器外形三視圖及網格圖見圖7.

圖7 Hyflex外形及網格圖Fig.7 Schematic and mesh of Hyflex

計算條件為:Ma=2.0，H=30 km，Sref=4.27 m2，Lref=4 m.計算攻角為:α =0°，5°，10°，15°，20°.減縮頻率定義為:k= ωLref/2u∞，本節計算時取 k=0.01，振幅取 αm=1°.

表2給出了各攻角下的直接阻尼導數、洗流時差導數及俯仰阻尼導數計算結果.可以看出，直接計算得到的俯仰阻尼導數與兩個分量相加得到的俯仰阻尼導數Cmq+最大相對誤差為6%，雖然沒有試驗數據作為對比，但通過兩種不同方法得到的俯仰阻尼導數相差不大，使得本文發展的阻尼導數計算方法的正確性得到了很好的自證.

表2 Hyflex阻尼導數計算結果Table2 Calculation results of pitch-damping derivatives for Hyflex

另外，從表2中可以看出，直接阻尼導數為負值，且絕對值隨攻角增加而單調增加;洗流時差導數也全為負值，且絕對值先增加后減小，隨攻角呈非線性變化，在－15°達到最大，在俯仰阻尼導數中所占比例不再是小量，反而起主導作用(最大比重達78%).與2.1節和2.2節研究的外形相比，Hyflex是帶翼升力體外形飛行器，尾部翼面的洗流時差效果更顯著.

3 結論

本文研究直接求解俯仰阻尼導數分量的方法，通過對HBS和基本帶翼導彈兩個標模外形及Hyflex飛行器進行數值模擬研究，結果表明:

1)采用的方法能夠準確預測飛行器外形的俯仰阻尼導數的各個分量，即使對于大攻角狀態也具有較好的預測精度，具備一定的工程實用價值;

2)對于彈丸類彈道物體，其在超聲速及高超聲速狀態下，洗流時差導數在俯仰阻尼導數中所占比例較小，但其符號可能大于零，起動不穩定作用，并且數值隨質心后移而增大;對于帶翼飛行器，超聲速狀態下，洗流時差導數在俯仰阻尼導數中所占比例較大.

直接阻尼導數不論從物理意義上來看還是數值計算都是負值，總起動穩定作用.

References)

［1］ 李周復.風洞特種試驗技術［M］.北京:航空工業出版社，2010:210-250.Li Z F.Special wind tunnel testing technology［M］.Beijing:Aviation Industry Press，2010:210-250(in Chinese).

［2］ Ren Y X.Evaluation of the stability derivatives using the sensitivity equations［J］.AIAA Journal，2008，46(4):912-917.

［3］ 劉偉，楊小亮，趙云飛.高超聲速飛行器加速度導數數值模擬［J］.空氣動力學學報，2010，28(4):426-429.Liu W，Yang X L，Zhao Y F.Numerical simulation of acceleration derivative of hypersonic aircraft［J］.Acta Aerodynamica Sinica，2010，28(4):426-429(in Chinese).

［4］ 孫海生.飛行器大攻角升沉平移加速度導數測量技術［J］.流體力學實驗與測量，2001，15(4):15-19.Sun H S.A measurement technique for derivatives of aircraft due to acceleration in heave and sideslip at high angle of attack［J］.Experiments and Measurements in Fluid Mechanics，2001，15(4):15-19(in Chinese).

［5］ Weinacht P.Navier-Stokes predictions of the individual components of the pitch-damping coefficient sum，ARL-TR-3169［R］.Adelphi:Army Research Laboratory，2004.

［6］ Weinacht P.Projectile performance，stability and free flight motion prediction using computational fluid dynamics［J］.Journal of Spacecraft and Rockets，2004，41(2):257-263.

［7］ Qin N，Ludlow D K，Shaw S T，et al.Calculation of pitch damping coefficients for projectiles，AIAA-1997-0405［R］.Reston:AIAA，1997.

［8］ 閻超.計算流體力學方法及應用［M］.北京:北京航空航天大學出版社，2006:18-25.Yan C.The methodology and application of computational fluid dynamics［M］.Beijing:Beihang University Press，2006:18-25(in Chinese).

［9］ McGowan G Z，Kurzen M J，Nance R P，et al.High fidelity approaches for pitch damping prediction at high angles of attack，AIAA-2012-2903［R］.Reston:AIAA，2012.

［10］ Hashimoto A，Hashizume M，Sunada S，et al.Unsteady analysis of aerodynamic derivatives on standard dynamics model，AIAA-2013-0343［R］.Reston:AIAA，2013.

［11］ East R A，Hutt G R.Comparison of predictions and experimental data for hypersonic pitching motion stability［J］.Journal of Spacecraft and Rockets，1988，25(3):225-233.

［12］ Uselton B L，Uselton J C.Test mechanism for measuring pitch damping derivatives of missile configurations at high angles of attack，AEDC-TR-75-43［R］.Tennessee:AEDC，1975.

［13］ Uselton B L，Jenke L M.Experimental missile pitch-and rolldamping characteristics at large angles of attack［J］.Journal of Spacecraft and Rockets，1977，14(4):241-247.

［14］ Murphy C H.Free flight motion of symmetric missiles，NO.1216［R］.Aberdeen:Army Ballistic Research Lab，1963.

［15］ Shirouzu M，Yamamoto M.Overview of the hyflex project，AIAA-1996-4524［R］.Reston:AIAA，1996.