基于CFD數(shù)值仿真技術(shù)的飛行器動導數(shù)計算

米百剛，詹浩，朱軍

（西北工業(yè)大學 翼型葉柵空氣動力國防科技重點實驗室，陜西 西安 710072）

基于CFD數(shù)值仿真技術(shù)的飛行器動導數(shù)計算

米百剛，詹浩，朱軍

（西北工業(yè)大學 翼型葉柵空氣動力國防科技重點實驗室，陜西 西安 710072）

為了研究表征飛行器動態(tài)穩(wěn)定特性的動導數(shù)，結(jié)合計算流體力學方法，基于滑移網(wǎng)格技術(shù)，采用小幅度強迫振動法和差分法兩種非定常方法，推導了俯仰組合動導數(shù)以及滾轉(zhuǎn)阻尼導數(shù)的計算表達式，并結(jié)合多參考系模型采用準定常方法計算滾轉(zhuǎn)阻尼導數(shù)。建立了基于計算流體力學技術(shù)的飛行器動導數(shù)計算方法。對國際動導數(shù)標模Einner導彈進行了計算驗證。兩種非定常方法對俯仰組合動導數(shù)的計算誤差分別為0．65%、6．13%，對滾轉(zhuǎn)阻尼導數(shù)的計算誤差分別為2．5%、0．06%。求解滾轉(zhuǎn)阻尼的準定常方法對于滾轉(zhuǎn)阻尼導數(shù)的計算誤差為2．67%。三種方法的計算結(jié)果與風洞試驗結(jié)果吻合的很好。非定常方法計算精度高，準定常方法可以快捷迅速得到結(jié)果。本文的數(shù)值方法不僅適用于超聲速，同樣可以用于亞、跨聲速時的動導數(shù)計算。

滑移網(wǎng)格；多參考系模型；動導數(shù)；計算流體力學技術(shù)

0 引 言

本文結(jié)合CED計算技術(shù)發(fā)展了飛行器超聲速動導數(shù)計算方法。以國外動導數(shù)標模［17］Einner導彈為例，基于滑移網(wǎng)格技術(shù)，對動態(tài)氣動設(shè)計特性中的滾轉(zhuǎn)阻尼導數(shù)，以及俯仰組合動導數(shù)進行了模擬，并采用多參考系方法，對滾轉(zhuǎn)阻尼導數(shù)進行了數(shù)值計算。計算結(jié)果表明，無論是基于滑移網(wǎng)格的非定常方法，還是結(jié)合多參考系模型的準定常方法，都具有比較高的計算精度。

1 計算方法

1.1 控制方程

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

其中V為任意控制體，W是守恒變量，F(xiàn)為無粘（對流）通矢量項，F(xiàn)v為粘性通量，?V為控制體的邊界，n為控制體邊界單位外法向矢量，Re為計算的雷諾數(shù)。

空間離散采用二階迎風格式——通量差分分裂（Roe-EDS）格式，采用雙時間步推進并用隱式LUSGS格式迭代求解，應(yīng)用當?shù)貢r間步長、殘值光順、預(yù)處理和多重網(wǎng)格加速收斂。湍流模型采用對逆壓梯度流動模擬精度較高的k-ωSST湍流模型。

1.2 滑移網(wǎng)格

滑移網(wǎng)格是一種特殊的動網(wǎng)格，它將問題的求解域劃分成多個計算域，即動區(qū)域和靜止區(qū)域。各計算域之間采用定義好的分界面進行關(guān)聯(lián)，動區(qū)域沿著滑移面進行運動。在滑移網(wǎng)格問題中，動區(qū)域運動是相對于靜止參考系進行跟蹤的，因此，沒有運動參考系附加在計算域上，簡化了穿過分界面的通量傳遞。與一般動網(wǎng)格相比，網(wǎng)格的運動僅沿著滑移面進行整體滑移，沒有網(wǎng)格變形重構(gòu)的過程，相對節(jié)省了產(chǎn)生新網(wǎng)格所需的計算量，而且運動時各區(qū)域網(wǎng)格質(zhì)量不發(fā)生變化，可以提高計算精度。

1.3 多參考系模型

多參考系模型（MRE）是不同旋轉(zhuǎn)或移動速度的每個區(qū)域的穩(wěn)態(tài)近似。同樣它將計算域分成多個小的子域，每個子域可以有自己的運動方式，或靜止或旋轉(zhuǎn)或平移。此方法適用于邊界上流動區(qū)域幾乎均勻混合的情況。以同時存在靜止和旋轉(zhuǎn)域情況為例，在兩個子域上分別建立靜止坐標系和旋轉(zhuǎn)坐標系，控制方程在每個子域分別求解，在交接面上通過將速度換算成絕對速度形式進行流場信息交換。由于進行的是穩(wěn)態(tài)近似，多參考系模型是一種很高效的方法。

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

1.4.1 強迫振動法

強迫飛行器做小幅度俯仰振蕩時，依據(jù)傅立葉公式展開，剛體飛行器所受非定常氣動力矩可以寫為表達式：

其中Mz是瞬態(tài)非定常氣動力矩值；Mz0是平衡位置處的氣動力矩是基頻諧波分量的氣動力矩幅值；ω是強迫飛行器振動頻率；λ是剛體飛行器強迫振動時位移和氣動力矩之間的相位差；u（t）為高次諧波分量。

依據(jù)氣動力導數(shù)的概念，上述小幅度強迫俯仰振動的剛體飛行器所受非定常氣動力矩還可以表示為：

當剛體飛行器以低頻ω做小幅度振蕩時，其模型運動方程可以簡化為：

展開式（2）并略去高次諧波分量，同時將式（4）代入略去高階分量，整理同類項可得：

當非定常問題計算足夠長時，令ωt＝2nπ，就可以抹去初始效應(yīng)的影響，氣動力矩達到一個周期性穩(wěn)態(tài)值，于是式（5）可以寫為：

式（6）即為基于傅立葉及泰勒展開方法推導出的俯仰組合動導數(shù)的表達式。

1.4.2 差分法

假定飛行器以兩種不同角速度ωz1，ωz2等速上仰運動至相同攻角，依據(jù)飛行力學小擾動理論對其力矩進行展開并略去高階量可得：

式（9）即為基于差分法得到的俯仰組合動導數(shù)公式。

1.4.3 準定常方法

前面提出的強迫振動法以及差分法都是非定常方法，可以用來求解縱向，橫航向動導數(shù)。而一般的飛行器，具有順流方向的面對稱性或者線對稱性，橫向氣動力在0°迎角，定速轉(zhuǎn)動狀態(tài)下一般保持不變，且橫向更多的是考慮穩(wěn)定滾轉(zhuǎn)特性，因而通常使用準定常方法進行滾轉(zhuǎn)動導數(shù)的計算。

飛行器進行小角速度等速滾轉(zhuǎn)時，令驅(qū)動力矩為Mx，阻尼力矩為Mkx，根據(jù)小擾動理論：

文中考慮零滾轉(zhuǎn)角下的動穩(wěn)定性導數(shù)，θ＝0，因此繞質(zhì)心轉(zhuǎn)動動力學方程簡化為：

注意到方程左端向為0，即可推導出滾轉(zhuǎn)動導數(shù)表達式：

其中Ix為飛行器對ox軸的轉(zhuǎn)動慣量。

1.5 邊界條件

本文使用的基于滑移網(wǎng)格以及多參考系模型的求解方法要求將流場劃分為動域以及靜域，兩個域之間通過交接面進行數(shù)據(jù)傳遞。動域中的網(wǎng)格剛化后沿著交接面進行滑移運動，以此來實現(xiàn)物面的運動。本文涉及到的邊界條件有遠場邊界，物面邊界以及交接面邊界。其中遠場邊界使用基于Riemann不變量的無反射邊界條件，物面邊界無滑移，絕熱，交接面使用基于通量守恒的交接面邊界。

2 動導數(shù)計算公式無量綱化

一般來講，計算動導數(shù)的公式都是無量綱的，因此，需要對本文的幾種方法進行無量綱化處理。

2.1 強迫振動法

根據(jù)國家標準，減縮頻率取為k＝ωl／2V*，其中l(wèi)為參考長度，V*為遠場自由來流速度，代入可得：

2.2 差分法

引入無量綱角速度，

帶入原式，

可得無量綱差分法公式：

2.3 準定常方法

3 動導數(shù)算例及分析

3.1 計算模型及網(wǎng)格

計算模型為國際標準動導數(shù)模型Einner導彈，圖1為計算用的Einner導彈模型及網(wǎng)格，其中d＝1m，采用結(jié)構(gòu)網(wǎng)格，根據(jù)本文的方法將流場分為動域和靜域，兩部分通過interface交接面連接，這個交接面的網(wǎng)格不要求完全一致，給復雜模型網(wǎng)格劃分帶來方便。最終劃分的網(wǎng)格總量為390萬，其中，動域網(wǎng)格量260萬，靜域130萬。

圖1 計算模型及網(wǎng)格Fig.1 Computational model and grid

3.2 算例驗證

3.2.1 強迫振動法求解俯仰組合動導數(shù)

使用非定常方法求解動導數(shù)的關(guān)鍵在于準確地預(yù)測瞬時力矩。本文選取M∞＝1．58，α0＝0°為計算狀態(tài)，減縮頻率為k＝0．0158226，振幅為αm＝1°，也即振蕩規(guī)律可寫為

圖2為計算得到的瞬時力矩隨迎角的變化曲線?？梢悦黠@地看出，本文的非定常強迫振動方法能夠很好地預(yù)測迎角與力矩系數(shù)之間的遲滯效應(yīng)，計算的俯仰力矩系數(shù)遲滯環(huán)與文獻［1］中的非常一致。

圖2 力矩系數(shù)隨迎角變化規(guī)律圖Fig.2 Time history of pitching moment coefficient with angle of attack in pitching cycle

根據(jù)前面得到的動導數(shù)無量綱計算公式，可求解俯仰組合動導數(shù)，計算結(jié)果見表1。

表1 俯仰組合動導數(shù)Table 1 Pitchingdynamic stabilityderivative by using small oscillation method

文獻［1］中的俯仰組合動導數(shù)計算值為-506，誤差為3．62%。由此可見本文的計算方法對俯仰組合動導數(shù)的計算是比較準確的。

一般使用強迫振動方法計算動導數(shù)多使用積分法進行結(jié)果辨識，本文則使用穩(wěn)定后的平衡位置處單點氣動數(shù)據(jù)計算動導數(shù)。表2是兩種方法的對比，可以看出，兩者均能比較準確的得到動導數(shù)結(jié)果，相比下，本文的單點數(shù)據(jù)計算量小，能夠方便快速得到結(jié)果。

表2 辨識方法對比Table 2 Identification methods of dynamic derivatives

為了進一步就本文的小幅度強迫振動方法進行驗證，分別計算在馬赫數(shù)為1．58，1．74，1．88，2．1，2．55五個狀態(tài)下的俯仰組合動導數(shù)。

圖3的結(jié)果表明，重心位置在Xg＝5．0d時，本文所計算的動導數(shù)與風洞試驗數(shù)據(jù)較為一致，反映出在超音速范圍內(nèi)，隨著馬赫數(shù)的增加，俯仰組合動導數(shù)的逐漸減小，且趨勢變緩和。由于目前的技術(shù)限制，動導數(shù)試驗結(jié)果的誤差帶相當大，與試驗結(jié)果相比，本文采用數(shù)值模擬技術(shù)較為準確地計算了Einner導彈的俯仰組合動導數(shù)，在文中給定的馬赫數(shù)范圍內(nèi)很好地模擬出動導數(shù)的變化趨勢。

圖3 不同馬赫數(shù)下的俯仰組合動導數(shù)Fig.3 Pitching dynamic stability derivative versus Mach number

3.2.2 差分法求解俯仰組合動導數(shù)

本文的計算初始迎角均為0°，導彈以5°／s，10°／s的角速度分別等速上仰至5°。然后根據(jù)得到的瞬時力矩即可求出俯仰組合動導數(shù)。計算結(jié)果見表3。

表3 差分法求解俯仰組合動導數(shù)Table 3 Pitching dynamic stability derivative by usingdifferential method

表3的結(jié)果表明，使用差分法也能比較準確地計算俯仰組合動導數(shù)。

3.2.3 準定常方法求解滾轉(zhuǎn)動導數(shù)

給定導彈以5°／s的速度繞軸滾轉(zhuǎn)，使用前面所述的準定常方法，計算得到的滾轉(zhuǎn)阻尼導數(shù)為：

同樣使用強迫振動法以及差分法求解滾轉(zhuǎn)動導數(shù)，并與準定常方法的結(jié)果進行比較，可得表4。

由表4可以看出，準定常雖然計算量小，耗時短，但是精度較低。而強迫振動以及差分方法精度比較高。對于動導數(shù)的數(shù)值計算，關(guān)鍵在于計算非定常氣動力，非定常方法基于這一點提出，因而精度也稍高，相對而言，準定常方法是一種近似的方法，對流場特性細節(jié)的捕捉能力稍弱，所以精度較低。

表4 滾轉(zhuǎn)動導數(shù)求解結(jié)果Table 4 Rolling damp dynamic derivative

4 結(jié) 論

采用基于滑移網(wǎng)格的非定常方法以及基于多參考系模型的準定常方法對Einner導彈標準模型動導數(shù)進行了計算研究和分析。相比以往的使用網(wǎng)格重構(gòu)進行計算的方法，滑移網(wǎng)格的引入使得計算量減小，并且由于不存在網(wǎng)格變形帶來的影響，因此可以提高求解的精度。由本文的算例可以看出所發(fā)展的非定常方法以及準定常方法能夠較為準確地計算俯仰及滾轉(zhuǎn)動導數(shù)。其中小幅度強迫振動法以及差分法兩種非定常方法的關(guān)鍵在于準確預(yù)測瞬時力矩；準定常方法雖然精度稍差，但是計算量小，快捷。同時，本文的三種方法不僅可以用于超音速動導數(shù)辨識，也同樣適用于亞音速以及跨音速的動導數(shù)計算。

本文采用的動導數(shù)辨識方法給出了與風洞試驗較為吻合的數(shù)值計算結(jié)果，以及一致的動導數(shù)的變化趨勢，表明所采用的動導數(shù)計算方法的正確性，具備較好的工程應(yīng)用前景。

［1］GREEN L L，SPENCE A M，MURPHY P C．Computation methods fordynamic stability and controlderivatives［R］．AIAA 2004-0015，2004．

［2］SHI A M．Numerical analysis of aircraft unsteady flow and aeroelastic problems［D］．Xi’an：Northwestern Polytechnical University，2005．（in Chinese）

史愛明．飛行器非定常流場和氣動彈性問題的數(shù)值分析研究［D］．西安：西北工業(yè)大學，2005．

［3］RONCH A D，VALLESPIN D．Computation ofdynamicderivatives using CED［R］．AIAA 2010-4817，2010．

［4］GREEN L L，SPENCE A M，MURPHY P C．Computation methods fordynamic stability and controlderivatives［R］．AIAA 2004-0015，2004．

［5］RONCH A D，GHOREYSHI M，BADCOCK K J，et al．Linear frequencydomain and harmonic balance predicitions ofdynamicderivatives［R］．AIAA 2010-4699，2010．

［6］PARK M A，GREEN L L．Steady-state computation of constant rotational ratedynamic stabilityderivatives［R］．AIAA 2000-4321，2000．

［7］KRAMER，BRIAN R．Experimental evaluation of superposition techniques applied todynamic aerodynamics［R］．AIAA 2002-0700，2002．

［8］ERDAL O，HSSAN U A．CED predicitions ofdynamicderivatives for missile［R］．AIAA 2002-0276，2002．

［9］SOO H P，YOONSIK K，JANG H K．Prediction ofdynamicdamping coefficients using unsteadydual-time stepping method［R］．AIAA 2002-0715，2002

［10］SCOTT M M，ELORET C．A reduced-frequency approach for calculatingdynamicderivatives［R］．AIAA 2005-840，2005．

［11］YUAN X X，ZHANG H X，XIE Y E．The pitching static／dynamicderivatives computation based on CED methods［J］．ACTA Aerod ynamica Sinica，2005，23（4）：458-463．（in Chinese）

袁先旭，張涵信，謝昱飛．基于CED方法的俯仰靜，動導數(shù)數(shù)值計算［J］．空氣動力學學報，2005，23（4）：458-463．

［12］SUN T，GAO Z H，HUANG J T．Identify of aircraftdynamicderivatives based on CED technology and analysis of reduce frequency［J］．Elight Dynamics，2011，29（4）：15-18．（in Chinese）

孫濤，高正紅，黃江濤．基于CED的動導數(shù)計算與減縮頻率影響分析［J］．飛行力學，2011，29（4）：15-18．

［13］SHI A M，YANG Y N，YE Z Y．A more accurate method for calculating transonicdynamicderivatives（TDDs）using present state-of-the-art CED［J］．Journal of Northwestern Polytechnical University，2008，26（1）：11-14．（in Chinese）

史愛明，楊永年，葉正寅．結(jié)合CED技術(shù)的跨音速動導數(shù)計算方法研究［J］．西北工業(yè)大學學報，2008，26（1）：11-14．

［14］LU X C，YE Z Y，ZHANG W W．A high efficient method for computingdynamicderivatives of supersonic／hypersonic aircraft［J］．Aeronautical Computing Technique，2008，38（3）：28-31．（in Chinese）

盧學成，葉正寅，張偉偉．超音速、高超聲速飛行器動導數(shù)的高效計算方法［J］．航空計算技術(shù)，2008，38（3）：28-31．

［15］SUN Z W，CHENG Z Y，BAI J Q，et al．A high efficient method for computingdynamicderivatives of aircraft based on quasisteady CED method［J］．Elight Dynamics，2010，28（2）：28-30．（in Chinese）

孫智偉，程澤蔭，白俊強，等．基于準定常的飛行器動導數(shù)的高效計算方法［J］．飛行力學，2010，28（2）：28-30．

［16］JIANG S J，LIU Y Q，DANG M L．A calculation method of aircraft roll-damping moment coefficientderivative based on steady NS equation［J］．Journal of Projectiles，Rockets，Missiles and Guidance，2008，28（1）：180-182．

將勝炬，劉玉琴，黨明利．基于定常NS方程的飛行器滾轉(zhuǎn)阻尼力矩系數(shù)導數(shù)計算方法［J］．彈箭與制導學報，2008，28（1）：180-182．

［17］GREENWELL D L．Erequency effects ondynamic stabilityderivatives obtained from small-amplitude oscillatory testing［J］．Journal of Aircraft，1998，35（5）：776-783．

Calculation ofdynamicderivatives for aircraft based on CFD technique

MI Baigang，ZHAN Hao，ZHU Jun
（National Key Laboratory of Aerodynamic Design and Research，Northwestern Polytechnical University，Xi’an 710072，China）

To study thedynamic stabilityderivatives of aircraft with CED technology，an unsteady method for calculating pitchingdynamic stabilityderivative and rollingdampdynamicderivative isdeveloped by using small amplitude oscillation model anddifferential model，and a quasi-steady method is also built to calculate rollingdampdynamicderivative；the unsteady method is based on sliding mesh technique while the quasi-steady method on both sliding mesh technique and moving reference frame model；compared to thedynamic mesh used before，the sliding mesh technique can improve the computational precision as it avoids the effect of meshdeformation．The finner missile is taken as example to calculate thedynamicderivatives．The errors of two unsteady methods to calculate longitudinal combinedderivatives are 0．65%and 6．13%respectively，while errors to lateral combinedderivatives are 2．5%and 0．06%，and the error of quasi-steady method to obtain lateral combinedderivative is 2．67%．The research shows the calculation result is consistent with the experimentaldata，although the quasi-steady method is less accurate than the unsteady one，nevertheless it can get result more quickly．Moreover，both the methods can be used for supersonic，subsonic and transonicdynamicderivative calculation．

sliding mesh；moving reference frame model；dynamicderivative；CED technology

V211.3

Adoi：10.7638／kqdlxxb-2012.0206

0258-1825（2014）06-0834-06

2012-12-13；

2013-02-25

米百剛（1989-），男，陜西富平人，博士研究生，研究方向：計算流體力學．E-mail:mibaigang＠m(xù)ail.nwpu.edu.cn

米百剛，詹浩，朱軍．基于CED數(shù)值仿真技術(shù)的飛行器動導數(shù)計算［J］．空氣動力學學報，2014，32（6）：834-839．

10．7638／kqdlxxb-2012．0206 MI B G，ZH AN H，ZHU J．Calculation ofdynamicderivatives for aircraft based on CED technique［J］．ACTA Aerodynamica Sinica，2014，32（6）：834-839．