黃宇， 閻超,*， 席柯， 王文

1.北京航空航天大學(xué) 航空科學(xué)與工程學(xué)院， 北京　100083 2.中國兵器工業(yè)導(dǎo)航與控制技術(shù)研究所, 北京　100089

?

基于數(shù)值虛擬飛行技術(shù)的飛行器動態(tài)特性分析

黃宇1， 閻超1,*， 席柯2， 王文1

1.北京航空航天大學(xué) 航空科學(xué)與工程學(xué)院， 北京100083 2.中國兵器工業(yè)導(dǎo)航與控制技術(shù)研究所, 北京100089

基于結(jié)合結(jié)構(gòu)重疊網(wǎng)格、閉環(huán)PID控制器、舵偏控制律、剛體六自由度運動和非定常Navier-Stokes方程求解等模塊的數(shù)值虛擬飛行技術(shù)，對“起源號”返回艙、基本帶翼導(dǎo)彈外形的多自由度非定常運動、受控特性及控制參數(shù)的整定開展了模擬。分析了不同自由度(DOF)下飛行器的運動特性，飛行器受擾動后的穩(wěn)定性及控制參數(shù)的整定。計算結(jié)果表明：利用數(shù)值虛擬飛行技術(shù)可有效地開展復(fù)雜流動下飛行物體非線性運動問題的研究，對研究飛行器在非線性流動下的動態(tài)特性、受控特性、流動機(jī)理研究以及控制律的設(shè)計檢驗具有較高的工程價值和應(yīng)用前景。

計算流體力學(xué)(CFD)； 重疊網(wǎng)格； PID控制； 虛擬飛行； 動導(dǎo)數(shù)； 控制參數(shù)整定； 動態(tài)特性； 非定常流動

飛行器在大迎角、過失速等狀態(tài)下飛行時，會出現(xiàn)復(fù)雜的分離渦等強非線性、非定常流動，與飛行器運動耦合將會使飛行器產(chǎn)生一系列的非定常運動現(xiàn)象，如極限環(huán)搖滾、俯仰振蕩等，從而引起飛行器氣動力的突變、遲滯及非對稱特性，極易導(dǎo)致飛行器運動進(jìn)入發(fā)散狀態(tài)，嚴(yán)重影響飛行器的安定性、穩(wěn)定性及受控性。由此，現(xiàn)代先進(jìn)空戰(zhàn)武器設(shè)計對飛行器的動態(tài)特性、穩(wěn)定性及控制效率等的預(yù)測技術(shù)提出了要求[1]。

飛行器出現(xiàn)非定常、非線性流動情況下，線性化數(shù)據(jù)不能正確地刻畫飛行器非線性的氣動力和運動規(guī)律[2-3]。為了能有效地研究飛行器非定常運動中的氣動、控制等問題，國內(nèi)外研究人員從20世紀(jì)80年代開始發(fā)展了一種新型的風(fēng)洞試驗技術(shù)，即風(fēng)洞虛擬飛行技術(shù)[4-6](Wind Tunnel Based Virtual Flight Testing)。通過將實驗縮比模型安裝有若干個自由度的支撐機(jī)構(gòu)上，在風(fēng)洞中開展運動模型的動力特性、控制特性的綜合實驗，在一定程度上模擬飛行器的真實飛行狀態(tài)，獲得飛行器動態(tài)品質(zhì)及控制特性。但由于風(fēng)洞試驗的固有缺陷，如支架/洞壁干擾、機(jī)械阻尼、相似準(zhǔn)則等難以完全滿足等情況，導(dǎo)致難以準(zhǔn)確地測量動態(tài)參數(shù)。此外，由于風(fēng)洞一般只能提供直勻流，且模型支撐機(jī)構(gòu)并不能在大范圍內(nèi)實現(xiàn)真實飛行自由度，因而模擬的仿真程度也受到制約。

近年來，計算流體力學(xué)(Computational Fluid Dynamics， CFD)及大規(guī)模高性能計算技術(shù)發(fā)展迅速，CFD在基礎(chǔ)理論和工程實踐上都取得了令人矚目的成就，以CFD技術(shù)為基礎(chǔ)實現(xiàn)的基于方程的數(shù)值虛擬飛行技術(shù)已經(jīng)成為可能[7-8]。數(shù)值虛擬飛行技術(shù)將飛行動力學(xué)方程/流體力學(xué)方程/控制律耦合實時求解，可開展全尺寸，真實運動/氣流參數(shù)的非定常模擬，避免了風(fēng)洞虛擬飛行技術(shù)的固有缺點，運動仿真程度更高。此外數(shù)值虛擬飛行技術(shù)靈活性更好、耗費低、周期短，在機(jī)理研究中對流場、運動參數(shù)的獲取更為方便和豐富，在實際工程應(yīng)用中可在各個階段介入研究，為設(shè)計及驗證飛行器氣動、控制等系統(tǒng)的特性提供有力的依據(jù)，填補風(fēng)洞虛擬飛行及飛行試驗的不足[9]。

本文基于數(shù)值虛擬飛行技術(shù)，開展了返回艙、帶翼導(dǎo)彈外形的動態(tài)特性研究。研究結(jié)果表明：數(shù)值虛擬飛行技術(shù)對復(fù)雜運動下的飛行器動態(tài)特性、飛行器的操縱特性、控制系統(tǒng)參數(shù)的整定等問題的研究具有較高的工程價值和應(yīng)用前景。

1　數(shù)值虛擬飛行方法

1.1流動主控方程

主控方程為非定常、可壓縮、雷諾平均的三維Navier-Stokes方程。采用來流參數(shù)為無量綱特征變量，對方程無量綱化，并將方程轉(zhuǎn)換到一般曲線坐標(biāo)系下，可得如下形式[8]：

(1)

式中：Q為守恒變量向量；F、G和H為對流通量；Fv、Gv和Hv為黏性通量；Re為來流雷諾數(shù)。

采用基于結(jié)構(gòu)網(wǎng)格的有限體積法求解，對流通量計算采用帶熵修正的Roe格式，采用MUSCL插值和vanAlbada限制器達(dá)到二階空間精度。時間推進(jìn)格式在定常計算時選用隱式LU-SGS方法，在非定常計算時選用雙時間步LU-SGS格式。采用一方程的Spalart-Allmaras湍流模型計算獲得湍流渦黏性。利用基于共享內(nèi)存的Open-MP并行技術(shù)加速計算。

1.2結(jié)構(gòu)重疊網(wǎng)格技術(shù)

網(wǎng)格生成是數(shù)值模擬的重要部分之一，好的網(wǎng)格質(zhì)量是獲得合理計算結(jié)果的前提條件。結(jié)構(gòu)重疊網(wǎng)格邏輯關(guān)系簡單，魯棒性好，適合開展工程計算。在大部分計算區(qū)域，結(jié)構(gòu)重疊網(wǎng)格計算方法同結(jié)構(gòu)網(wǎng)格完全一樣，并能保證和結(jié)構(gòu)網(wǎng)格相同的計算精度、計算效率及流動特性捕捉能力，而網(wǎng)格間的重疊技術(shù)降低了復(fù)雜拓?fù)渫庑蔚木W(wǎng)格生成難度[10-12]。另一方面，虛擬模擬實際問題時，被模擬物體子部件間往往會出現(xiàn)相對運動，如舵面/襟翼等部件受控偏轉(zhuǎn)、飛行器脫離掛架等，結(jié)構(gòu)重疊網(wǎng)格將各子部件劃為獨立的子區(qū)域，可有效地處理這些多體運動、分離問題[12]。求解多體、部件運動等位置或形狀隨時間變化的問題時，網(wǎng)格無需重新生成，各子區(qū)域網(wǎng)格隨部件作剛性運動，每一時間步只需要進(jìn)行一次重疊過程以建立插值信息關(guān)系，相比網(wǎng)格重構(gòu)、變形等動網(wǎng)格方法實現(xiàn)更為簡單高效，且適用于任意幅度的運動[12-13]。

1.3飛行器剛體6DOF運動方程組

飛行器或運動部件作為剛體處理,在空間中剛體運動有6個自由度(DOF)，剛體在空間的平動、轉(zhuǎn)動、位置及姿態(tài)的變化可以用6個動力學(xué)方程和6個運動學(xué)方程充分描述[14]。

剛體質(zhì)心動力學(xué)方程為

(2)

式中:a和V分別為慣性系下加速度和速度矢量；P和F分別為推力和氣動力在體軸坐標(biāo)系下矢量，下標(biāo)b表示體軸系；G為慣性系下重力矢量；L為坐標(biāo)轉(zhuǎn)換矩陣；m為質(zhì)量。

慣性坐標(biāo)系下質(zhì)心運動學(xué)方程為

(3)

體軸坐標(biāo)系下剛體繞質(zhì)心轉(zhuǎn)動的動力學(xué)方程為

(4)

(5)

由式(4)可得體軸系下繞質(zhì)心轉(zhuǎn)動的動力學(xué)方程為

(6)

再結(jié)合體軸系下角速度分量同姿態(tài)角變化關(guān)系式(7)，可獲得關(guān)于俯仰角θ、偏航角φ及滾轉(zhuǎn)角γ的剛體動力學(xué)和運動學(xué)聯(lián)立的方程組為

(7)

1.4控制系統(tǒng)?；?/p>

飛行控制系統(tǒng)一般由舵回路、穩(wěn)定回路和制導(dǎo)回路這3個反饋回路構(gòu)成。這3個回路主要實現(xiàn)氣動力(氣動面)控制、擔(dān)負(fù)穩(wěn)定與控制功能，以及按照一定導(dǎo)引律給出控制指令[15]。

在數(shù)值虛擬飛行中可以使用任意的控制器，復(fù)雜的控制器的實現(xiàn)、考察不在本文研究范圍中。為簡化起見，本文采用PID控制器，實現(xiàn)舵偏反饋控制，并形成穩(wěn)定的制導(dǎo)回路來跟蹤指定迎角變化規(guī)律。PID控制器如圖1所示。

圖1　PID控制器框圖Fig.1　Block diagram of PID controller

圖中：R(t)為設(shè)定值，C(t)為檢測值，E(t)為設(shè)定值和檢測值間的偏差，可見PID控制器是通過對偏差進(jìn)行運算獲得控制值U(t)。PID運算包含3部分，第1部分為比例控制Proportional，比例控制器為一個增益(放大)器，其將輸入信號放大Kp倍，Kp稱為增益(放大)系數(shù)，增加Kp可提高系統(tǒng)準(zhǔn)確性，降低系統(tǒng)慣性，但會影響控制穩(wěn)定性。第2部分為積分控制器Integral，增加積分控制器參數(shù)Ki可以減小控制系統(tǒng)的穩(wěn)態(tài)誤差，但有可能引起系統(tǒng)振蕩。第3部分為比例微分控制器Derivative，改變Kd可影響系統(tǒng)響應(yīng)速度，即影響控制系統(tǒng)的快速性[16]。

在模擬中，考慮到非定常CFD計算的物理時間步長遠(yuǎn)遠(yuǎn)小于舵機(jī)的響應(yīng)時間，因此假設(shè)舵面每隔Td接受一次控制信號。在Td時間內(nèi)，舵偏規(guī)律設(shè)為

(8)

式中：Δd為需要接受控制指令后舵面需要偏轉(zhuǎn)的總舵偏角；d(t)為舵偏角關(guān)于時間的函數(shù)；int為取整函數(shù)。

1.5氣動/運動/控制集成流程

數(shù)值虛擬飛行的核心流程是氣動/運動/控制(CFD/RBD/FCS)耦合求解。其中，CFD模塊采用非定常方法求解獲得每一物理時間步飛行器受到的氣動力和力矩。剛體運動(RigidBodyDynamics，RBD)模塊根據(jù)力和力矩情況，求解剛體6DOF方程，更新物體位置(網(wǎng)格坐標(biāo))。

本文使用松耦合算法[17]求解剛體動力學(xué)方程與流動控制方程的耦合問題，采用龍格庫塔法求解動力學(xué)方程，采用雙歐法求解飛行器姿態(tài)變化。在獲得剛體下一時刻的位移和姿態(tài)角后，由式(9)獲得n+1時刻的網(wǎng)格點坐標(biāo)為

Rn+1=Ln+1·R0

(9)

式中：Ln+1為由n+1時刻的姿態(tài)角確定的轉(zhuǎn)換矩陣；R0為初始的網(wǎng)格坐標(biāo)。

控制系統(tǒng)(FlightControlSystem，F(xiàn)CS)模塊根據(jù)飛行器姿態(tài)與期望的偏差形成反饋控制，結(jié)合重疊網(wǎng)格等動網(wǎng)格技術(shù)控制氣動操縱面，控制飛行器姿態(tài)。然后再反復(fù)上述過程，直到運動結(jié)束。流程如下：

1) 為減小計算時間，采用定常計算獲得初始流場。

2) 以初場為基礎(chǔ)開展非定常流動計算。

3) 每推進(jìn)一個物理時間步，求解飛行器的氣動力和力矩，代入RBD求解獲得物體位移，更新整體網(wǎng)格位置。

4) 根據(jù)FCS，控制操縱面運動，更新操縱面物體(網(wǎng)格)位置。

5) 開展下一物理時間步流場計算，重復(fù)步驟2～5，直至物理時間推進(jìn)結(jié)束。

2　數(shù)值方法驗證

用三維的超聲速彈道外形[18](HyperBallisticShape，HBS)對程序和方法的正確性和可信度進(jìn)行驗證。HBS外形如圖2所示，具有試驗數(shù)據(jù)可供對比，是一個廣泛用于超聲速動態(tài)非定常計算程序驗證的標(biāo)準(zhǔn)算例。

圖2　超聲速彈道模型外形及尺寸Fig.2　Shape and size of HBS model

采用文獻(xiàn)[18]中來流參數(shù)開展計算，具體為：來流馬赫數(shù)6.85，質(zhì)心在彈長0.72倍處，以前段直徑d為特征長度的雷諾數(shù)為0.72×106。對比在0°、5°、10°、15°和20° 迎角下計算和試驗獲得的俯仰導(dǎo)數(shù)，俯仰動導(dǎo)數(shù)采用振幅為1° 的正弦強迫運動計算獲得，強迫運動的減縮頻率取0.01，計算物理時間步長為2ms。圖3(a)為計算獲得的靜導(dǎo)數(shù)Cmα同試驗值對比，程序準(zhǔn)確計算出了約在8° 迎角附近導(dǎo)數(shù)符號的改變。圖3(b)為俯仰動導(dǎo)數(shù)同試驗值對比，計算同試驗數(shù)據(jù)符合較好，表明本文程序和非定常方法是可行的。

圖3　HBS俯仰靜、動導(dǎo)數(shù)計算與試驗數(shù)據(jù)對比Fig.3　Comparison of computational static and dynamic pitch derivative with test data of HBS

3　數(shù)值模擬結(jié)果及分析

3.1返回艙外形平面運動動態(tài)特性研究

起源號返回艙網(wǎng)格如圖4所示，具體外形尺寸見文獻(xiàn)[19]，模型物理參數(shù)見表1。

圖4　起源號返回艙表面及對稱面網(wǎng)格Fig.4　Surface and symmetry plane grid of Genesis capsule

表1　起源號返回艙模型物理屬性

Table 1　Physical property of Genesis capsule model

ItemValued/mm69.42M/g633Izz/kg·m23.01×10-4Iyy/kg·m21.99×10-4L/mm44.33Xcg/mm23.83

圖5為初始馬赫數(shù)Ma0為4.5，釋放迎角為0.5° 時返回艙總迎角隨射程變化規(guī)律同文獻(xiàn)數(shù)據(jù)對比。由圖可見，該狀態(tài)彈道的前段總迎角較小，在120 m后迎角開始發(fā)散，計算也出現(xiàn)了相同的趨勢。需要說明的是，文獻(xiàn)中并未說明返回艙釋放的角速度，計算時設(shè)初始角速度為零，從結(jié)果來看這不會對后段彈道迎角發(fā)散這一規(guī)律產(chǎn)生較大影響。圖6為Ma=3.5，大氣參數(shù)為海平面時，采用CFD自由振蕩計算獲得的動導(dǎo)數(shù)同文獻(xiàn)[19]中數(shù)據(jù)對比。計算結(jié)果同文獻(xiàn)數(shù)據(jù)接近，符合該返回艙在馬赫數(shù)大于2.5時在大迎角下為動穩(wěn)定的結(jié)論。

圖5　Ma0=4.5，釋放迎角為0.5° 時，返回艙總迎角隨彈道變化規(guī)律計算與試驗對比Fig.5　Comparison of computational trajectory history of total angle of attack with test data at Ma0=4.5 and release α=0.5°

圖6　Ma=3.5時，返回艙俯仰動導(dǎo)數(shù)計算同文獻(xiàn)對比Fig.6　Comparison of computational capsule dynamic pitch derivative with reference data at Ma=3.5

以上計算結(jié)果同試驗結(jié)果吻合良好，可確定本文方法的正確性。下面的研究主要關(guān)注不同自由度運動對返回艙運動特性的影響，因此本外形算例來流參數(shù)統(tǒng)一為：來流馬赫數(shù)為3.5，飛行高度為25 km，來流壓力、溫度根據(jù)地球大氣參數(shù)表獲得，計算皆采用層流假設(shè)，物理時間步長統(tǒng)一取為0.3 ms。由于自由俯仰運動時，該系統(tǒng)是自治系統(tǒng)，其極限環(huán)狀態(tài)(振幅和頻率)與初值無關(guān)，因此，為了較快出現(xiàn)極限環(huán)運動，計算預(yù)置10° 迎角，并且設(shè)置返回艙的質(zhì)量為表1中的1/4，以減小慣量，使運動迅速出現(xiàn)振蕩，以節(jié)約計算時間。

3.1.1平面單自由度俯仰運動

當(dāng)返回艙繞質(zhì)心做單自由度俯仰振動時，其動力學(xué)方程為

(10)

在某迎角的小的鄰近范圍內(nèi)，可以假設(shè)俯仰靜、動導(dǎo)數(shù)為常量，則式(10)為一典型的二階阻尼系統(tǒng)，其解析解為

(11)

式中：ξl為阻尼；ωl為振動頻率，具體表達(dá)式為

(12)

圖7為返回艙的靜、動導(dǎo)數(shù)隨迎角變化情況。由圖可見，各迎角下該外形靜導(dǎo)數(shù)變化不大，而動導(dǎo)數(shù)在小迎角處是正值，表現(xiàn)為動不穩(wěn)定；迎角大于5° 后，動導(dǎo)數(shù)變?yōu)樨?fù)值，表現(xiàn)為動穩(wěn)定。因此，返回艙在小迎角下受到擾動，由于動不穩(wěn)定性，擾動振幅將增大，在大迎角時，氣流起阻尼作用，從而衰減振幅。在一定條件下，這樣的動導(dǎo)數(shù)變化規(guī)律可使得飛行器出現(xiàn)俯仰極限環(huán)振蕩狀態(tài)。

圖7　“起源號”返回艙各迎角下俯仰靜、動導(dǎo)數(shù)Fig.7　Static and dynamic pitch derivative of Genesiscapsule under different angles of attack

圖8為自由俯仰振蕩迎角隨時間的變化規(guī)律?？梢姡汕笆鰟訉?dǎo)數(shù)特性，從10° 迎角釋放外形后，飛行器振幅逐漸增加，當(dāng)振幅達(dá)到15° 左右時，飛行器出現(xiàn)等周期的極限環(huán)振蕩。

圖8　自由俯仰運動迎角隨時間變化規(guī)律Fig.8　Time history of angle of attack of free pitching motion

3.1.2平面俯仰+沉浮雙自由度運動

下面研究引入質(zhì)心沉浮運動后，返回艙的運動特性。圖 9為雙自由度和單自由度運動下飛行器迎角隨時間變化規(guī)律。由圖可見，引入沉浮運動后，返回艙仍形成等周期和振幅極限環(huán)俯仰振動，但相比單自由度振動，其振幅加大，頻率減小。

圖9　雙自由度和單自由度運動下迎角隨時間變化規(guī)律的對比Fig.9　Comparison of time histories of angle of attack of one and two DOF motion

圖10為返回艙質(zhì)心沉浮速度V隨時間變化規(guī)律，其與迎角隨時間變化有1/4周期的相位差。此外，注意到，由沉浮運動引起的迎角變化不會超過0.02°，這表明引起俯仰振動頻率和振幅變化的主要原因并不是沉浮速度，而是引入沉浮運動后，系統(tǒng)阻尼變小導(dǎo)致的，這跟文獻(xiàn)中理論分析一致[20-21]。返回艙平面雙自由度運動下的俯仰靜、動導(dǎo)數(shù)及失穩(wěn)特性和單自由度運動的類似，在此不再贅述。

圖10　返回艙雙自由度運動質(zhì)心沉浮速度隨時間的變化規(guī)律Fig.10　Time history of centroid velocity of two DOF motion of return capsule

3.1.3平面自由飛三自由度運動

返回艙返回時，在氣動阻力的作用下減速運動，引入流向運動自由度，即可考察返回艙平面自由飛行的動態(tài)特性，這種自由飛更加接近真實飛行情況。圖11為平面自由飛運動下飛行器迎角隨時間的變化規(guī)律。由圖可見，考慮氣動阻力后，返回艙振幅更大，這表明飛行器變得更為動不穩(wěn)定。此外，可以看出，振蕩頻率逐漸減小，振蕩周期變長，這是由于氣動阻力使飛行器速度降低造成的。

圖11　平面自由飛和雙自由度下迎角隨時間的變化規(guī)律對比Fig.11　Comparison of time histories of angle of attack of planar free and two DOF motion

圖12為平面自由飛運動下返回艙質(zhì)心沉浮速度隨時間的變化規(guī)律?？梢?，自由飛下沉浮速度振幅比雙自由度運動下的振幅有所增加，且周期逐漸變長。圖13為質(zhì)心流向速度隨時間的變化規(guī)律，空氣阻力的存在使飛行器流向速度由3.5馬赫降到了2.2馬赫左右。同樣，基于自由飛計算開展參數(shù)辨識，也可獲得相應(yīng)的靜、動導(dǎo)數(shù)，在此不再累述。

圖12　平面自由飛質(zhì)心沉浮速度隨時間變化規(guī)律Fig.12　Time history of centroid velocity of planar free motion

圖13　平面自由飛質(zhì)心流向速度隨時間變化規(guī)律Fig.13　Time history of flow direction centroid velocity of planar free motion

3.2受控導(dǎo)彈縱向數(shù)值虛擬飛行模擬

計算外形為尾部安裝有4個翼的導(dǎo)彈。該外形具有大量的靜、動態(tài)試驗數(shù)據(jù)及計算對比數(shù)據(jù)[22-23]。該導(dǎo)彈外形涵蓋了大部分戰(zhàn)術(shù)導(dǎo)彈的典型氣動布局特征。導(dǎo)彈為十字翼布局，尾舵為均勻厚度。計算外形尺寸及舵號定義見圖14(圖中CG為模型重心)，模型物理參數(shù)見表 2。圖15為Ma=4.42，Red=2.96×107時采用不同時間步長計算獲得的俯仰角隨射程變化規(guī)律同試驗值之間的對比，其中Red是以彈身直徑d為特征長度的雷諾數(shù)；圖 16為Red=8.6×106，Ma=1.96，無量綱質(zhì)心位置X/d=6.1時導(dǎo)彈各迎角下動導(dǎo)數(shù)計算與試驗對比。以上計算結(jié)果同試驗結(jié)果吻合良好，可確定本文方法的正確性。此外由圖 15可見，采用0.5 ms和1 ms時間步長的計算結(jié)果同試驗值得差異皆不大，因此為減少計算機(jī)時，之后計算中物理時間步長均取1 ms。本文只考慮俯仰通道模擬，因此只考慮1號和2號舵實現(xiàn)舵偏功能，舵后緣向下偏轉(zhuǎn)定義為正舵偏，采用重疊網(wǎng)格實現(xiàn)舵偏功能。此外，在下述計算皆采用一方程的Spalart-Allmaras湍流模型模擬湍流黏性，子迭代收斂殘差取為0.01。

圖14　基本帶翼導(dǎo)彈模型外形及尺寸Fig.14　Shape and size of basic finner

表2　基本帶翼導(dǎo)彈計算模型物理屬性

Table 2Physical property of calculation model of basic finner

ItemValueL/m0.3M/kg1.58Ixx/kg·m21.92×10-4Iyy/kg·m297.85×10-4d/m0.03Xcg/m0.165

圖15　導(dǎo)彈俯仰角隨射程變化計算同試驗值對比Fig.15　Comparison of computational projectile pitch angle variation with range and test data

圖16　Red=8.6×106時動導(dǎo)數(shù)計算試驗值對比Fig.16　Comparison of computational dynamic derivative with test data at Red=8.6×106

3.2.1飛行器姿態(tài)穩(wěn)定特性研究

飛行器的穩(wěn)定特性包括動穩(wěn)定性和動操縱性。當(dāng)飛行器受外界擾動后的運動體現(xiàn)其動穩(wěn)定性，當(dāng)飛行器受操縱力控制后的運動體現(xiàn)其動操縱性。設(shè)計擾動規(guī)律，開展飛行器姿態(tài)穩(wěn)定研究，可以檢驗飛行器不同的運動模態(tài)，測試飛行器在各種極端情況下的靜、動穩(wěn)定性，且可為氣動參數(shù)辨識提供足夠的信息。

采用時域法分析一個系統(tǒng)的性能，一般采用一些典型的輸入作為測試信號，如階躍輸入，脈沖輸入等[16]，結(jié)合導(dǎo)彈飛行的特點，本文選取如下4種舵偏規(guī)律作為擾動輸入，測試基本外形導(dǎo)彈在不同舵偏擾動下的姿態(tài)穩(wěn)定性：

1) 階躍型舵偏輸入

d(t)=5°t≥0

(13)

2) 方形脈沖舵偏輸入

(14)

3) 梯形舵偏輸入

(15)

4) 正弦舵偏輸入

(16)

式中:d(t)為舵偏角隨時間的函數(shù);t為物理時間。下面計算中來流條件為Ma=2，高度為10 km，來流壓力及密度根據(jù)大氣參數(shù)表獲得，飛行器初始迎角為0°,無姿態(tài)角速度及姿態(tài)角加速度。

基于小擾動理論，基本外形導(dǎo)彈帶舵面控制的無量綱縱向運動方程可寫為

(17)

式中：Iz為轉(zhuǎn)動慣量；θ為俯仰角。

式(17)的特征方程為

(18)

圖17為階躍舵偏下導(dǎo)彈俯仰角響應(yīng)曲線及俯仰角速度響應(yīng)曲線。階躍舵偏反映了舵偏(擾動)突變的一種極限情況，而飛行受擾動后收斂于平衡位置的過程反映了飛行的動穩(wěn)定性。由圖可見，階躍舵偏(擾動)引起的角速度達(dá)200°/s左右，飛行器俯仰角振蕩收斂到-5.5° 左右，振蕩次數(shù)、幅度衰減情況反映了飛行器的操縱品質(zhì)。式(17)在階躍舵偏下的解為

(19)

式中：αeq為配平攻角；φ為相位裕度，具體表達(dá)式為

(20)

結(jié)合式(18)與式(19)，即可辨識出俯仰阻尼導(dǎo)數(shù)等參數(shù)。

圖18為偶極脈沖舵偏下導(dǎo)彈俯仰角響應(yīng)曲線及俯仰角速度響應(yīng)曲線。在t=0.2 s時，此時導(dǎo)彈具有正俯仰角速度，0.2 s時刻后，舵偏變負(fù)，產(chǎn)生的抬頭力矩加劇了抬頭的趨勢。當(dāng)干擾頻率與飛行器固有頻率接近時，可能引起飛行器共振。偶極脈沖舵偏可反映飛行器受到瞬時、隨機(jī)且頻率接近飛行器固有頻率的干擾的動態(tài)性能。

圖19為梯形脈沖舵偏下導(dǎo)彈俯仰角響應(yīng)曲線及俯仰角速度響應(yīng)曲線。梯形脈沖舵偏反映了等速舵偏與階躍舵偏的組合情況。由圖可見，此時超調(diào)量小于階躍型舵偏，且振動頻率也小于階躍型舵偏，最大俯仰角速度約為100°/s，這種舵偏規(guī)律代表了緩慢的機(jī)動飛行。

圖17　階躍舵偏下的俯仰響應(yīng)曲線Fig.17　Pitch response curve of phase step impulse rudder deflexion

圖18　偶極脈沖舵偏下的俯仰響應(yīng)曲線Fig.18　Pitch response curve of dipole impulse rudder deflexion

圖19　梯形脈沖舵偏下的俯仰響應(yīng)曲線Fig.19　Pitch response curve of echelon impulse rudder deflexion

圖20為諧波脈沖舵偏下導(dǎo)彈俯仰角響應(yīng)曲線及俯仰角速度響應(yīng)曲線。由圖可見，此種情況下，舵偏(擾動)引起的最大俯仰角速度及俯仰力矩相對前3種情況是最小的，此時飛行器實施精確跟蹤，處于理想化的反復(fù)修正過程。

圖20　諧波脈沖舵偏下的俯仰響應(yīng)曲線Fig.20　Pitch response curve of harmonic wave impulse rudder deflexion

3.2.2基于數(shù)值虛擬飛行對控制參數(shù)的整定

飛行器在飛行過程中，控制系統(tǒng)會根據(jù)飛行狀態(tài)發(fā)出實時控制指令，修正飛行姿態(tài)。系統(tǒng)受控后，要么導(dǎo)致受控發(fā)散，要么趨于受控狀態(tài)。飛行器是否會出現(xiàn)受控發(fā)散，即控制是否具有穩(wěn)定性，可以采用3.3.1節(jié)中的方法檢驗。而對于受控穩(wěn)定情況，則要考慮控制系統(tǒng)的準(zhǔn)確性，快速性及平穩(wěn)性。

控制系統(tǒng)的準(zhǔn)確性，快速性及平穩(wěn)性與控制系統(tǒng)的參數(shù)有很強的關(guān)系。同一種控制方式，面對各種各樣的飛行器，滿足上述要求的控制參數(shù)也大不相同。即使對于同一個飛行器，由于飛行過程中燃料消耗導(dǎo)致的質(zhì)量、慣性矩及重心的物理參數(shù)的改變，也會使飛行器受控后響應(yīng)發(fā)生改變。而對于處于非定常/非線性流動下的飛行器，氣動力、力矩的突變、遲滯特性更會使得基于線化氣動導(dǎo)數(shù)等參數(shù)獲得控制參數(shù)不能在全范圍內(nèi)保證控制的穩(wěn)定、準(zhǔn)確、快速及平穩(wěn)性。因此，整定控制參數(shù)，是數(shù)值虛擬飛行一個重要的應(yīng)用，也是風(fēng)洞等常規(guī)手段難以實現(xiàn)的。

本文采用PID控制器，實現(xiàn)舵偏反饋控制，跟蹤指定迎角規(guī)律。測試了不同PID參數(shù)(如表3 所示)下導(dǎo)彈的縱向受控性能，以獲得較好的控制參數(shù)。

模擬的來流條件為：Ma=2，高度為10 km，初始迎角為a=0°，初始俯仰角速度為零，角加速度也為零。導(dǎo)彈姿態(tài)角變化要求為：機(jī)動到5° 迎角飛行，超調(diào)不超過20%，并維持姿態(tài)穩(wěn)定。按上述要求及計算模型物理屬性，根據(jù)PID控制理論[16]，選取如下范圍的PID控制參數(shù)及舵受控時間Td，對該范圍的參數(shù)整定，以達(dá)到最佳控制效果。

表3　PID控制參數(shù)

圖21～圖24為各控制參數(shù)下導(dǎo)彈的響應(yīng)時間歷程。從圖中可以看到，控制系統(tǒng)根據(jù)飛行器姿態(tài)角的誤差發(fā)出控制指令，控制舵面，飛行器實施機(jī)動，然后反復(fù)調(diào)節(jié)，直至飛行器達(dá)到期望姿態(tài)。由相圖可見，初始階段，飛行器的俯仰角速度增加，迎角也隨之增加，當(dāng)迎角接近指定值5° 時，微分控制器開始制動，減小俯仰角速度，當(dāng)姿態(tài)角度超過5° 后，在控制器作用下，舵面反復(fù)調(diào)整，使得俯仰角速度變負(fù)，最后逐漸收斂于相圖上的螺旋點。

對比Case 1和Case 2可以看出，Case 1超調(diào)量約20%，Case 2超調(diào)量約10%，Case 1的最大俯仰角速率約為33°/s，而Case 2為30°/s。這表明增加Kp，會提高系統(tǒng)的響應(yīng)速度，但會使得超調(diào)增加，振蕩過渡過程變長。對比Case 3和Case 4相圖可以看出，在迎角達(dá)到2° 左右時，控制系統(tǒng)開始制動，即減小俯仰角速率。制動后Case 3的最大俯仰角速率為39°/s，Case 4的最大俯仰角速率為42°/s。這表明增加Kd，使得系統(tǒng)的制動提前，降低俯仰角速度峰值，能減少過渡振蕩過程。對比Case 2和Case 3，Td減小加快了俯仰角速率，提高了系統(tǒng)響應(yīng)速度。

圖21　Case 1的俯仰響應(yīng)曲線Fig.21　Pitch response curve of Case 1

圖22　Case 2的俯仰響應(yīng)曲線Fig.22　Pitch response curve of Case 2

圖23　Case 3的俯仰響應(yīng)曲線Fig.23　Pitch response curve of Case 3

上述算例表明，Case 3的控制參數(shù)有較好的特性。下面采用Case 3控制參數(shù)對有初始俯仰角速度的狀態(tài)進(jìn)行控制，驗證控制效果，初始俯仰角速度設(shè)為50°/s。

圖25為采用Case 3的控制參數(shù)下的俯仰角控制曲線和相圖。導(dǎo)彈在初始俯仰角速度的作用下，迎角增加。由于導(dǎo)彈是靜、動穩(wěn)定，且此時舵偏角較小，俯仰角速度迅速減小變?yōu)樨?fù)值。當(dāng)舵偏角負(fù)向變大時，產(chǎn)生抬頭力矩，迎角增加，然后逐漸調(diào)整到預(yù)設(shè)迎角。初始俯仰角速度的存在使得彈體動阻尼效應(yīng)較為明顯，且初始角速度越大，動態(tài)效應(yīng)表現(xiàn)越強。

圖24　Case 4的俯仰響應(yīng)曲線Fig.24　Pitch response curve of Case 4

圖25　有初始俯仰角速度時的俯仰響應(yīng)曲線Fig.25　Pitch response curve of state with initial pitching angular speed

4　結(jié)　論

本文利用數(shù)值虛擬飛行技術(shù)，對返回艙外形及受控基本帶翼導(dǎo)彈外形開展了動態(tài)特性研究。

對返回艙外形的研究表明：返回艙外形具有靜態(tài)穩(wěn)定性和動不穩(wěn)定性。對于平面三自由度俯仰運動，俯仰阻尼導(dǎo)數(shù)決定了返回艙的動穩(wěn)定性。引入沉浮運動使得返回艙動穩(wěn)定性變差。引入流向自由度后，氣動阻尼使飛行器速度降低，動穩(wěn)定性變差，且極限環(huán)周期逐漸變長。

對受控基本外形導(dǎo)彈的研究表明：通過輸入典型的擾動(舵偏)信號，可以檢測導(dǎo)彈系統(tǒng)的穩(wěn)定、受控性能。不同控制參數(shù)對系統(tǒng)的控制性能有較大影響。增大增益系數(shù)，可以提高響應(yīng)速度，但會增加超調(diào)量，引起振蕩。增加微分系數(shù)，可以控制系統(tǒng)制動幅度及時間，減小超調(diào)量，使振蕩過渡過程盡快結(jié)束。

上述研究表明，結(jié)合結(jié)構(gòu)重疊網(wǎng)格技術(shù)、閉環(huán)控制律的數(shù)值虛擬飛行技術(shù)能夠有效地研究常規(guī)手段無法研究的飛行器動態(tài)特性問題。其應(yīng)用靈活性好，周期短，對探索飛行器氣動/運動耦合問題，控制參數(shù)的整定，控制律的研究和檢測，獲取飛行器非定常的運動導(dǎo)數(shù)以及飛行的安全性測試具有較好的工程應(yīng)用價值及應(yīng)用前景。

[1]呂光男. 風(fēng)洞虛擬飛行實驗中的飛行力學(xué)與控制問題研究[D]. 南京: 南京航空航天大學(xué), 2009： 1-13

LYU G N. The problem of flight mechanics and control in wind tunnel based virtual flight testing[D]. Nanjing: Nanjing University of Aeronautics and Astonautics, 2009:1-13 (in Chinese).

[2]李周復(fù). 風(fēng)洞特種實驗技術(shù)[M]. 北京: 航空工業(yè)出版社, 2010: 210-250.

LI Z F. Special wind tunnel testing technology[M]. Beijing: Aviation Industry Press, 2010: 210-250 (in Chinese).

[3]劉偉, 楊小亮, 趙云飛. 高超聲速飛行器加速度導(dǎo)數(shù)數(shù)值模擬[J]. 空氣動力學(xué)學(xué)報, 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]RATLIFF C L, MARQUART E J. An assessment of a potential test technique: virtual flight testing (VFT): AIAA-1995-3415[R]. Reston: AIAA, 1995.

[5]GEBERT G, KELLY J, LOPEZ J. Wind tunnel based virtual flight testing: AIAA-2000-0829[R], Reston: AIAA, 2000.

[6]MAGILL J C, CATALDI P, MORENCY J R, et al. Demonstration of a wire suspension for wind-tunnel virtual flight testing[J]. Journal of Spacecraft and Rockets, 2009, 46(3): 624-633.

[7]張來平, 馬戎, 常興華, 等. 虛擬飛行中氣動運動和控制耦合的數(shù)值模擬技術(shù)[J]. 力學(xué)進(jìn)展, 2014, 44: 201410.

ZHANG L P, MA R, CHANG X H, et al. Review of aerodynamics/kinematics/flight-control coupling methods in virtual flight simulations[J]. Advances in Mechanics, 2014, 44: 201410.

[8]閻超. 計算流體力學(xué)方法及應(yīng)用[M]. 北京: 北京航空航天大學(xué)出版社, 2006： 18-24.

YAN C. Computational fluid dynamics methods and applications. Beijing[M]. Beijing: Beihang University Press 2006: 18-24 (in Chinese).

[9]SALAS M D. Digital flight: The last CFD aeronautical grand challenge[J]. Journal of Scientific Computing, 2006, 28(2-3): 479-505.

[10]ROGERS S E,ROTH K,NASH S M,et al.Advances in overset CFD processes applied to subsonic high-lift aircraft: AIAA-2000-4216[R]. Reston: AIAA, 2000.

[11]GOETZ H K, JAMES E K, HENRY C L, et al. Validation of overflow for computing plume effects during the ARES I stage separation process: AIAA-2011-0170[R]. Reston: AIAA, 2011

[12]范晶晶, 閻超, 張輝. 重疊網(wǎng)格洞面優(yōu)化技術(shù)的改進(jìn)與應(yīng)用[J]. 航空學(xué)報, 2010, 31(6): 1127-1133.

FAN J J, YAN C, ZHANG H. Improvement of hole-surface optimization technique in overset grids and its application[J]. Acta Aeronauticaet Astronautica sinica, 2010, 31(6)： 1127-1133 (in Chinese).

[13]袁武. 新型重疊網(wǎng)格方法研究及在復(fù)雜多體氣動問題中的應(yīng)用[D]. 北京： 北京航空航天大學(xué), 2012: 24-35.

YUAN W. Investigations on novel chimera grid methods and its applications to complex multibody aerodynamic problems[D]. Beijing: Beihang University, 2012: 24-35 (in Chinese).

[14]肖業(yè)倫， 航空航天器運動的建模—飛行動力學(xué)的理論基礎(chǔ)[M]. 北京： 北京航空航天大學(xué)出版社， 2003: 64-78.

XIAO Y L. Modeling of the motion of space vehicle—the basic of flight mechnics[M]. Beijing: Beihang University Press 2003: 64-78 (in Chinese).

[15]吳森堂, 費玉華. 飛行控制系統(tǒng)[M]. 北京: 北京航空航天大學(xué)出版社, 2009: 85-93.

WU S T, FEI Y H, Flight control system[M]. Beijing: Beihang University Press, 2009: 85-93 (in Chinese).

[16]MEI X R, ZHUANG X Y. Principle of auto control[M]. Beijing: Science Press, 2009: 54-61.

[17]栗長江. 高超聲速飛行器多自由度動態(tài)特性研究[D]. 長沙: 國防科學(xué)技術(shù)大學(xué), 2010: 51-54.

LI C J. Study of mutidimensional motion of hypersonic flight vehicle[D]. Changsha: National University of Science and Technology, 2010: 51-54 (in Chinese).

[18]EAST R A, HUTT G R. Comparison of predictions and experimental data for hypersonic pitching motion stability[J]. Journal of Spacecraft, 1988, 25(3): 225-233

[19]CHEATWOOD F, WINCHENBACH G, HATHAWAY W, et al. Dynamic stability testing of the genesis sample return capsule: AIAA-2000-1009[R]. Reston: AIAA, 2000.

[20]袁先旭. 非定常流動數(shù)值模擬及飛行器動態(tài)特性分析研究[D]. 綿陽: 中國空氣動力研究與發(fā)展中心研究生部, 2002: 101-117.

YUAN X X. Numerical simulation for unsteady flowsand research on dynamic characteristicsof vehicle[D]. Mianyang: China Aerodynamics Research and Development Center. 2002: 101-117 (in Chinese).

[21]席柯. 虛擬飛行數(shù)值模擬方法及飛行器動態(tài)特性研究[D]. 北京： 北京航空航天大學(xué), 2015.

XI K. Virtual flight numerical simulation method and research on dynamic characteristics of vehicle[D]. Beijing, Beihang University, 2015 (in Chinese).

[22]DUPUIS A D, HATHAWAY W. Aeroballistic range tests of the basic finner reference projectile at supersonic velocities: DREV-TM-9703[R]. Valcartier: Defense Research Establishment, 1997.

[23]BHAGWANDIN V A, SAHU J. Numerical prediction of pitch damping stability derivatives for finned projectiles: AIAA-2011-3028[R]. Reston: AIAA, 2011.

黃宇男， 博士研究生。主要研究方向： 計算流體力學(xué)、 混合網(wǎng)格數(shù)值方法。

Tel: 010-82338071

E-mail: huangyu@buaa.edu.cn

閻超男， 博士， 教授， 博士生導(dǎo)師。主要研究方向： 計算流體力學(xué)、 空氣動力學(xué)。

Tel: 010-82317019

E-mail: yanchao@buaa.edu.cn

Analysis of flying vehicle’s dynamic characteristics based onnumerical virtual flight technology

HUANG Yu1, YAN Chao1,*, XI Ke2, WANG Wen1

1. School of Aeronautic Science and Engineering, Beihang University, Beijing100083, China 2. Institute of Ordnance Industry Navigation and Control Technology, Beijing100089, China

With the numerical virtual flight technology based on structured overlapping grid, closed PID controller, rudder control, rigid body motion and unsteady N-S equation solver, the unsteady motion characteristic, dynamic stability and adjustment of control parameter of Genesis capsule and basic finner projectile have been simulated. The motion characteristic with different degree of freedom (DOF) motion, the motion stability and the adjustment of control parameter with disturbance of those flying vehicles have been analyzed. The calculation results show that the study of flying vehicle’s nonlinear motion under the condition of complex fluid flows can be effectively carried out by numerical virtual flight technology, which has practical value and application prospect in the area of simulation and prediction of flying vehicle’s motion and controlling characteristic, fluid mechanics study and design of control rule under the condition of unsteady and nonlinear aerodynamics.

CFD; overlapping grid; PID control; virtual flight; dynamic derivative; control parameters adjustment; dynamic characteristic; unsteady flow

2016-01-11； Revised： 2016-02-17； Accepted： 2016-04-06； Published online： 2016-04-2015:26

. Tel.： 010-82317019E-mail: yanchao@buaa.edu.cn

2016-01-11； 退修日期： 2016-02-17； 錄用日期： 2016-04-06；

時間： 2016-04-2015：26

www.cnki.net/kcms/detail/11.1929.V.20160420.1526.004.html

.Tel.： 010-82317019E-mail:yanchao@buaa.edu.cn

10.7527/S1000-6893.2016.0120

V211.3

A

1000-6893(2016)08-2525-14

引用格式： 黃宇， 閻超， 席柯， 等． 基于數(shù)值虛擬飛行技術(shù)的飛行器動態(tài)特性分析[J]． 航空學(xué)報, 2016, 37(8): 2525-2538. HUANG Y, YAN C, XI K, et al. Analysis of flying vehicle’s dynamic characteristics based on numerical virtual flight technology[J]. Acta Aeronautica et Astronautica Sinica, 2016, 37(8): 2525-2538.

http://hkxb.buaa.edu.cnhkxb@buaa.edu.cn

URL： www.cnki.net/kcms/detail/11.1929.V.20160420.1526.004.html