達(dá)興亞，周為群，趙忠良，陶 洋

（中國空氣動力研究與發(fā)展中心 高速空氣動力研究所，綿陽 622661）

0 引 言

高機(jī)動性導(dǎo)彈突出強(qiáng)調(diào)大迎角快速轉(zhuǎn)彎能力，極小展弦比窄條翼在大迎角不存在渦破裂現(xiàn)象，縱向氣動特性極為優(yōu)異，是高機(jī)動性導(dǎo)彈經(jīng)常采用的布局形式［1］（圖1）。但是，窄條翼導(dǎo)彈強(qiáng)烈的非線性會帶來嚴(yán)重的翼／舵氣動干擾，導(dǎo)致全彈氣動特性隨滾轉(zhuǎn)角和迎角急劇變化，在機(jī)動時極易誘發(fā)非指令的自激搖滾振蕩。

圖1 窄條翼示意圖（NASA和德國IRIS－T）Fig.1 Strake missiles（NASA and German IRIS－T）

搖滾是飛行器滾轉(zhuǎn)方向的自由振蕩，由于以滾轉(zhuǎn)方向?yàn)橹?，一般簡化為單自由度運(yùn)動來開展研究［2］。通常有三種途徑：一是通過風(fēng)洞實(shí)驗(yàn)獲取搖滾特性；二是通過數(shù)值模擬搖滾運(yùn)動，分析流場變化，從微觀上搞清搖滾的流動機(jī)制；三是通過非線性動力學(xué)理論，分析搖滾運(yùn)動的動力學(xué)特性，搞清搖滾的產(chǎn)生機(jī)理，給出失穩(wěn)判據(jù)和預(yù)測方法。張涵信等人指出［3］，搞清搖滾的產(chǎn)生機(jī)理是十分必要的：即什么條件下?lián)u滾是穩(wěn)定的，什么條件下演化為周期搖滾。由于國內(nèi)研究起步晚，對這類導(dǎo)彈的搖滾機(jī)理還了解不多。

在過去的工作中［4］，對窄條翼布局導(dǎo)彈模型進(jìn)行了風(fēng)洞搖滾實(shí)驗(yàn)，發(fā)現(xiàn)窄條翼導(dǎo)彈模型具有兩個搖滾區(qū)間，第一搖滾區(qū)間在20°迎角附近，迎角范圍極?。坏诙u滾區(qū)間在35°迎角，迎角范圍較寬。實(shí)驗(yàn)同時明確了搖滾運(yùn)動與窄條翼／尾舵干擾密切相關(guān)，但尚不清楚對應(yīng)的動力學(xué)特性。本文針對模型第二搖滾區(qū)開展數(shù)值計算，重點(diǎn)分析搖滾的動力學(xué)特性，搞清搖滾的產(chǎn)生機(jī)理。

1 數(shù)值計算方法

1.1 控制方程

搖滾數(shù)值模擬常使用的控制方程有速度勢方程［5］、歐拉方程［6］，以及 NS方程［7］。對三角翼研究表明，速度勢方程和歐拉方程能成功模擬出極限環(huán)振蕩，同時也表明搖滾可能主要是非定常無粘現(xiàn)象主導(dǎo)［8］。三維NS方程對分離預(yù)測更準(zhǔn)確，能得到與實(shí)驗(yàn)更接近的結(jié)果。本文流場控制方程采用三維非定常NS方程：

式中Q為守恒變量，F(xiàn)、G、H為無粘矢通量，F(xiàn)v、Gv、Hv為粘性矢通量。轉(zhuǎn)動方程為定義在彈體坐標(biāo)系的剛體運(yùn)動方程：

其中φ為滾轉(zhuǎn)角，Ix為轉(zhuǎn)動慣量，ω為角速度，M為力矩，上標(biāo)b表示彈體坐標(biāo)系下的分量。

1.2 計算方法［9］

流場求解器是基于結(jié)構(gòu)網(wǎng)格的有限體積求解器，粘性項(xiàng)采用James中心差分，無黏項(xiàng)采用Roe格式，以F為例：

各參數(shù)具體意義參見文獻(xiàn)［10］。時間推進(jìn)采用雙時間步法［11］，湍流模型選取S－A 模型［12］。

轉(zhuǎn)動方程使用Adams預(yù)估校正法求解，氣動／運(yùn)動耦合計算采用雙時間步Adams預(yù)估校正策略，具體思路是：在CFD完成n時刻計算并得到氣動力后，采用Adams－Bashforth預(yù)測格式預(yù)測n＋1時刻的滾轉(zhuǎn)角和滾轉(zhuǎn)角速度，更新姿態(tài)和CFD邊界條件；之后，CFD采用雙時間步法開始計算n＋1時刻流場，在雙時間步的每一內(nèi)迭代步結(jié)束后，CFD得到新的氣動力，此時采用Adams－Moulton格式修正滾轉(zhuǎn)角和角速度，在內(nèi)迭代結(jié)束或者狀態(tài)變量收斂后，推進(jìn)到n＋2時刻。

在雙時間步內(nèi)迭代過程實(shí)現(xiàn)隱式修正，三階Adams－Moulton格式為：這里n代表時間步，N代表內(nèi)迭代步，如果計算收斂，隨著內(nèi)迭代步的推進(jìn)，xn＋1，N＋1將趨于常值。

研究表明，以上三階耦合計算方法在保證流場收斂和一定計算精度條件下，可以顯著增大時間步長，縮短計算時間，詳細(xì)分析參見文獻(xiàn)［9］。

2 搖滾運(yùn)動數(shù)值模擬結(jié)果及討論

2.1 計算網(wǎng)格及狀態(tài)

計算模型為鈍頭／窄條翼／尾舵組合導(dǎo)彈模型（類似IRIS－T，如圖1），網(wǎng)格為標(biāo)準(zhǔn)多塊結(jié)構(gòu)對接形式，在模型的周向分布了223個網(wǎng)格點(diǎn)，尾舵網(wǎng)格進(jìn)行了加密處理，網(wǎng)格總量約580萬，其上游距頭部10L，下游距后緣10L，遠(yuǎn)場邊界距中心線7L，L為全彈長度。為提高并行計算效率，將網(wǎng)格分為111塊，最大塊網(wǎng)格量約15萬，保證每個計算核心分配到大致相等的計算量。

計算針對窄條翼導(dǎo)彈模型的第二搖滾區(qū)，Ma＝0.6，α＝35°，“×”字型狀態(tài)，Re和Ix取風(fēng)洞試驗(yàn)對應(yīng)值。計算采取預(yù)偏滾轉(zhuǎn)角φ0、初始角速度為0的起動策略。

圖2 計算網(wǎng)格空間拓?fù)銯ig.2 Space topo of computational grid

2.2 計算結(jié)果及討論

由圖3可見，模型從φ0＝5°自由釋放，在經(jīng)歷短暫的幾個振蕩周期后，模型很快進(jìn)入極限環(huán)振蕩（Limit Cycle Oscillation，LCO），振幅約10.14°，頻率20Hz。圖4力矩遲滯曲線形狀上比較“肥胖”，說明遲滯現(xiàn)象明顯，模型吸收和耗散能量的速度很快，導(dǎo)致模型很快進(jìn)入極限環(huán)狀態(tài)。圖5相位圖可以看出，進(jìn)入極限環(huán)振蕩后，新周期的相位與前一周期完全重合。

圖3 滾轉(zhuǎn)角時間歷程Fig.3 Time history of roll angle

圖4 遲滯曲線Fig.4 Hysteresis loop

圖5 相平面Fig.5 Phase plane

此狀態(tài)下對應(yīng)的風(fēng)洞試驗(yàn)結(jié)果為：振動均方差σγ＝8°，頻率14.3Hz，計算結(jié)果與試驗(yàn)結(jié)果振幅吻合較好，而頻率比試驗(yàn)值大5.7Hz。引起頻率差異主要有兩方面原因：一是在氣動載荷的作用下，模型轉(zhuǎn)動時可能存在較大的機(jī)械阻尼，目前在計算中還沒有一種普適的方法引入機(jī)械阻尼；二是模型轉(zhuǎn)動時連軸承一起轉(zhuǎn)動，計算時沒有計入這部分轉(zhuǎn)動慣量。

極限環(huán)顯著的特征是一個周期內(nèi)吸收與耗散的能量相等，總能量變化為0［13］。當(dāng)模型處于低能狀態(tài)時，會不斷吸收能量而趨向極限環(huán)，當(dāng)處于高能狀態(tài)時，會不斷耗散能量而趨向極限環(huán)。圖5中由內(nèi)向外發(fā)展的相位圖即是由低能態(tài)向極限環(huán)過渡的過程，圖4中的內(nèi)環(huán)為順時針方向，系統(tǒng)吸收能量，兩邊環(huán)是逆時針方向，能量被耗散，其面積和等于順時針環(huán)面積，將吸收的能量耗散掉，一個周期內(nèi)系統(tǒng)能量變化為0，振蕩既不發(fā)散，也不收斂，維持極限環(huán)狀態(tài)。

3 搖滾動力學(xué)特性分析

3.1 搖滾產(chǎn)生的動力學(xué)機(jī)理

圖6給出了迎風(fēng)面和背風(fēng)面窄條翼和尾舵的遲滯曲線。從斜率上看，背風(fēng)舵和背風(fēng)翼的斜率為負(fù)，且斜率較大，迎風(fēng)舵遲滯環(huán)雖大，但斜率小，迎風(fēng)翼的斜率為正，因此，導(dǎo)彈的橫向靜穩(wěn)定性主要由背風(fēng)舵和背風(fēng)翼提供。

四條曲線中僅迎風(fēng)舵力矩為順時針環(huán)，即有能量注入，運(yùn)動發(fā)散；背風(fēng)翼和背風(fēng)舵為逆時針環(huán)，系統(tǒng)對外做功，能量被耗散。迎風(fēng)翼也為逆時針單環(huán)，雖然面積較小，但對系統(tǒng)的動穩(wěn)定性貢獻(xiàn)為正。也可以看出，迎風(fēng)舵在平衡滾轉(zhuǎn)角處的遲滯環(huán)面積很大，而背風(fēng)舵的遲滯環(huán)面積很小，導(dǎo)致背風(fēng)舵不能抵消迎風(fēng)舵的動不穩(wěn)定性，使模型喪失滾轉(zhuǎn)阻尼，吸收能量，運(yùn)動逐漸發(fā)散；隨著振幅的增大，背風(fēng)舵的動穩(wěn)定性增強(qiáng)，而迎風(fēng)舵的動不穩(wěn)定性降低，整個模型從動不穩(wěn)定變?yōu)閯臃€(wěn)定，耗散能量，運(yùn)動受到抑制，最終進(jìn)行等幅等周期極限環(huán)振蕩。

圖6 窄條翼和尾舵遲滯曲線Fig.6 Hysteresis loops for strakes and fins

3.2 初始滾轉(zhuǎn)角對搖滾的影響

大多數(shù)情況下，極限環(huán)是穩(wěn)定的，模型從任意非平衡狀態(tài)都將穩(wěn)定到該極限環(huán)［14］；特別情況下極限環(huán)是不穩(wěn)定的，系統(tǒng)還存在其它極限環(huán)，模型在擾動作用下會在極限環(huán)之間來回跳動，進(jìn)入混沌運(yùn)動［15］。

取初始滾轉(zhuǎn)角γ0＝15°，模型同樣進(jìn)入極限環(huán)振蕩，如圖7所示，相圖由外向內(nèi)發(fā)展，即由高能態(tài)過渡到極限環(huán)。該極限環(huán)與γ0＝5°極限環(huán)基本重合，說明確實(shí)存在一個穩(wěn)定的極限環(huán)，模型在高／低能狀態(tài)或初始微擾動作用下，都將進(jìn)入該極限環(huán)振蕩。在模擬時未考慮外部擾動，如果外部擾動過大，極限環(huán)的穩(wěn)定性不再起作用，也可能進(jìn)入混沌運(yùn)動，這種情況不在本文討論范圍。

圖7 初始滾轉(zhuǎn)角的影響Fig.7 Effect of initial roll angle

3.3 轉(zhuǎn)動慣量對搖滾的影響

搖滾實(shí)驗(yàn)采用了鋁制模型，未作轉(zhuǎn)動慣量影響研究，為了明確轉(zhuǎn)動慣量的影響，假設(shè)采用鋼制模型（密度約為鋁的3倍），將轉(zhuǎn)動慣量增大至3倍后進(jìn)行模擬。模擬的振幅9.3°，頻率9.4Hz，振幅變化不大，頻率變化明顯。

由圖8中力矩遲滯曲線可見，力矩遲滯曲線都呈雙“8”環(huán)，環(huán)的大小基本一致，只是兩邊的小環(huán)稍有差異。比較圖5和圖9，增大轉(zhuǎn)動慣量導(dǎo)致最大滾轉(zhuǎn)角速度從1200°／s降低到600°／s，說明在這一較寬廣的轉(zhuǎn)速范圍，氣動和運(yùn)動的耦合特性基本相似，角速度對力遲滯特性和流場跟隨性的影響較小。

圖8 轉(zhuǎn)動慣量的影響Fig.8 Effect of moment of inertia

圖9 相平面（Ix＝3I0）Fig.9 Phase plane（Ix＝3I0）

4 結(jié) 論

通過對窄條翼導(dǎo)彈模型的搖滾運(yùn)動模擬，得到了搖滾的動力學(xué)特性，結(jié)論如下：

1）多個模擬狀態(tài)表明，窄條翼導(dǎo)彈模型力矩遲滯曲線呈雙“8”環(huán)，中間環(huán)為順時針方向，兩邊環(huán)為逆時針方向；

2）窄條翼導(dǎo)彈模型搖滾的動力學(xué)原因是：由于背風(fēng)舵不能提供足夠的動穩(wěn)定性，模型在迎風(fēng)舵的動不穩(wěn)定作用下運(yùn)動發(fā)散，隨著振幅的增大，背風(fēng)舵動穩(wěn)定性增強(qiáng)，迎風(fēng)舵動不穩(wěn)定性減弱，從而抑制了振蕩發(fā)散，最終進(jìn)入極限環(huán)振蕩；

3）第二搖滾區(qū)的極限環(huán)是穩(wěn)定的，模型在任意初始狀態(tài)或微擾動作用下都將進(jìn)入該極限環(huán)振蕩；

4）在非定常特性較強(qiáng)時，存在一個轉(zhuǎn)動慣量范圍，搖滾振幅變化不大，振蕩周期變化明顯。

［1］ALLEN J M.Aerodynamics of an axisymmetric missile concept having cruciform strakes and in－line tail fins from mach 0.60to 4.63［R］.NASA TM 2005－213541，2005.

［2］KATZ J.Wing／vortex interactions and wing rock［J］.Progress inAerospaceSciences，1999，35：727－750.

［3］ZHANG H X，LIU W，XIE Y F，et al.On the rocking motion and its dynamic evolution of a swept delta wing［J］.ACTA AerodynamicaSinica，2006，24（1）：5－9.（in Chinese）張涵信，劉偉，謝昱飛，等.后掠三角翼的搖滾及其動態(tài)演化問題［J］.空氣動力學(xué)學(xué)報，2006，24（1）：5－9.

［4］DA X Y，ZHAO Z L，TAO Y，et al.Experimental investigation on free－to－roll motion of strake missile experimental［J］.JournalofExperimentsinFluidMechanics，2012，26（6）：40－44.（in Chinese）達(dá)興亞，趙忠良，陶洋，等.窄條翼導(dǎo)彈模型搖滾特性試驗(yàn)研究［J］.實(shí)驗(yàn)流體力學(xué)，2012，26（6）：40－44.

［5］KONSTADINOPOULOS P，MOOK D T，NAYFEH A H.Subsonic wing rock of slender delta wings［J］.JournalofAircraft，1985，22（3）：223－228.

［6］KANDIL O A，SALMAN A A.Three－dimensionla simulation of slender delta wing rock and divergence［R］.AIAA－92－0280，1992.

［7］LIU W，ZHANG H，ZHAO H.Numerical simulation and physical characteristics analysis for slender wing rock［J］.JournalofAircraft，2006，43（3）：858－861.

［8］ARENA A，NELSON R.A discrete vortex model for prediction wing rock of slender wings［R］.AIAA－92－4497，1992

［9］DA X Y，TAO Y，ZHAO Z L.Research on virtual flight simulation based on prediction－correction coupling method and chimera grid［J］.ActaAeronauticaetAstronauticaSinica，2012，33（6）：977－983 （in Chinese）達(dá)興亞，陶洋，趙忠良.基于預(yù)估校正和嵌套網(wǎng)格的虛擬飛行數(shù)值模擬［J］.航空學(xué)報，2012，33（6）：977－983.

［10］ROE P.Approximate Riemann solvers，parameter vectors，and difference schemes［J］.JournalofComputationalPhysics，1981，43（2）：357－372.

［11］DUBUC L，et al.Solution of the unsteady Euler equations using an implicit dual－time method［J］.AIAAJournal，1998，36（8）：1417－1424.

［12］SPALART P R，ALLMARAS S.A one－equation turbulence model for aerodynamic flows［R］.AIAA－92－0439，1992.

［13］NELSON R C，PELLETIER A.The unsteady aerodynamics of slender wings and aircraft undergoing large amplitude maneuvers［J］.ProgressinAerospaceSciences，2003，39：185－248.

［14］LIU Z L.Control and simulation investigation of wing rock phenomenon［D］.［PhD thesis］.Concordia University，Canada，2004.

［15］RONG Z，DENG X Y，WANG B，et al.Experimental investigation on flow characteristics during double limit cycle oscillation induced by forebody vortex［J］.JournalofExperimentsinFluid Mechanics，2010，24（3）：19－24.（in Chinese）榮臻，鄧學(xué)鎣，王兵，等.前體渦誘導(dǎo)雙極限環(huán)搖滾流動特性的實(shí)驗(yàn)研究［J］.實(shí)驗(yàn)流體力學(xué)，2010，24（3）：19－24.