楊文，尚祖銘，卜忱

(中航工業(yè)空氣動(dòng)力研究院氣動(dòng)發(fā)展部，黑龍江哈爾濱150001)

0 引言

隨著戰(zhàn)斗機(jī)飛行包線擴(kuò)大到過失速區(qū)域及各種復(fù)雜機(jī)動(dòng)的完成，作用在飛機(jī)上的氣動(dòng)力呈現(xiàn)高度的非線性及時(shí)間相關(guān)性，常規(guī)的準(zhǔn)定常氣動(dòng)力模型已不再適用。建立準(zhǔn)確、實(shí)用的大迎角非定常氣動(dòng)力模型是進(jìn)行現(xiàn)代飛機(jī)飛行動(dòng)力學(xué)分析、飛控系統(tǒng)設(shè)計(jì)及飛行仿真研究的前題。

最早的氣動(dòng)力模型是Bryan提出的線性模型。然而在大迎角下，該模型很難描述非線性氣動(dòng)力特性。Tobak［1］采用階躍響應(yīng)的方法建立了氣動(dòng)力系數(shù)與飛行狀態(tài)參數(shù)變化歷程的積分關(guān)系式，給出了氣動(dòng)力系數(shù)非線性、非定常數(shù)學(xué)模型的一般表達(dá)。雖然這種模型的理論較為完備，但過于復(fù)雜，實(shí)際應(yīng)用中很難使用。Goman等［2］首先提出了狀態(tài)空間模型的概念，它引用表征流動(dòng)分離點(diǎn)或渦破裂位置的內(nèi)部狀態(tài)變量，將氣動(dòng)力系數(shù)與流動(dòng)現(xiàn)象聯(lián)系起來，不僅表達(dá)式簡單，易于應(yīng)用實(shí)現(xiàn)，而且有清晰的物理意義，便于對大迎角非定常氣動(dòng)力特性進(jìn)行分析。但這種模型最大的缺陷在于橫航向機(jī)動(dòng)過程中還很難實(shí)際應(yīng)用。近年來，一種所謂的“黑箱”方法也引起了許多學(xué)者的興趣，它避開了復(fù)雜的物理機(jī)理，直接從飛行試驗(yàn)或者風(fēng)洞試驗(yàn)中獲得氣動(dòng)力數(shù)據(jù)，結(jié)合相應(yīng)的運(yùn)動(dòng)狀態(tài)變量，采用新型數(shù)學(xué)研究方法建立氣動(dòng)力與飛行狀態(tài)之間的關(guān)系式。史志偉等［3］應(yīng)用模糊邏輯方法建立了飛機(jī)非定常氣動(dòng)力模型，但這種純數(shù)學(xué)方法在建模樣本外的預(yù)測能力不太理想。

目前國內(nèi)外關(guān)于狀態(tài)空間方法的公開文獻(xiàn)主要是針對縱向運(yùn)動(dòng)建立模型，而對橫航向機(jī)動(dòng)的研究很少涉及。本文基于某戰(zhàn)斗機(jī)的大幅滾轉(zhuǎn)振蕩試驗(yàn)數(shù)據(jù)，應(yīng)用升力面分解方法［4］探索狀態(tài)空間模型在橫向機(jī)動(dòng)過程中的推廣應(yīng)用，并利用試驗(yàn)數(shù)據(jù)驗(yàn)證了所建立的非定常氣動(dòng)力模型。

1 試驗(yàn)內(nèi)容

1.1 試驗(yàn)系統(tǒng)

試驗(yàn)是在中航工業(yè)空氣動(dòng)力研究院FL-8回流式閉口低速風(fēng)洞中進(jìn)行的。風(fēng)洞試驗(yàn)段長5.5 m，截面尺寸3.5 m×2.5 m，最大風(fēng)速73 m/s。

動(dòng)態(tài)試驗(yàn)平臺(tái)采用新設(shè)計(jì)改造的雙自由度大幅振蕩試驗(yàn)系統(tǒng)，圖1為雙自由度振蕩系統(tǒng)結(jié)構(gòu)圖。

圖1 雙自由度振蕩系統(tǒng)Fig.1 Two-degree-of-freedom oscillation system

模型正裝時(shí)液壓馬達(dá)驅(qū)動(dòng)彎刀使模型進(jìn)行俯仰運(yùn)動(dòng)，伺服電機(jī)能夠驅(qū)動(dòng)模型滾轉(zhuǎn)振蕩。運(yùn)動(dòng)過程中，模型的實(shí)際迎角和側(cè)滑角可以通過機(jī)構(gòu)角轉(zhuǎn)換得到:

式中:?為平臺(tái)轉(zhuǎn)盤角;υ為大彎刀轉(zhuǎn)角;φ為尾撐支桿滾轉(zhuǎn)角。圖2為模型采用尾撐的某飛機(jī)大幅振蕩試驗(yàn)照片。

圖2 動(dòng)態(tài)試驗(yàn)照片F(xiàn)ig.2 Dynamic test photo

1.2 測控系統(tǒng)及數(shù)據(jù)處理

試驗(yàn)采用動(dòng)態(tài)專用六分量內(nèi)式天平，天平測得的電信號(hào)由FL-8風(fēng)洞VXI采集系統(tǒng)采集，數(shù)字信號(hào)由工控機(jī)處理，試驗(yàn)原始數(shù)據(jù)由大幅振蕩試驗(yàn)處理程序處理成體系下的體軸六分量系數(shù)。動(dòng)態(tài)濾波采用傅里葉變換數(shù)字濾波器，試驗(yàn)現(xiàn)場可以實(shí)時(shí)顯示曲線。

2 非定常橫向狀態(tài)空間模型

2.1 狀態(tài)空間模型

對于現(xiàn)代大多數(shù)飛行器而言，翼型都被設(shè)計(jì)成湍流邊界層。湍流邊界層的一個(gè)特點(diǎn)是流動(dòng)分離一般從后緣開始并隨著翼型迎角的增大而向前緣移動(dòng)。在這種流動(dòng)類型中，設(shè)有一個(gè)無量綱坐標(biāo)x∈［0，1］，x=0對應(yīng)于完全分離流場，x=1對應(yīng)于完全附著流。參數(shù)x就是Goman提出的與流動(dòng)分離相關(guān)的狀態(tài)變量［2］，對于小展弦比的戰(zhàn)斗機(jī)，它表征機(jī)翼上表面渦破裂的位置變量，并且滿足如下的常微分方程:

式中:x0(αeff)為狀態(tài)變量關(guān)于有效迎角的函數(shù)，且文獻(xiàn)［5］提出了一個(gè)補(bǔ)充變量y，且y=1－x，新的定義使得狀態(tài)變量的物理意義更加清晰。在新的定義下，完全分離流或完全渦破裂對應(yīng)于y=1，這個(gè)補(bǔ)充的狀態(tài)變量雖然不能描述翼型上分離流或機(jī)翼上渦破裂發(fā)生的具體位置，但能夠描述分離流或渦破裂引起的效應(yīng)，如果將上述微分方程中的x用y來代替，則:

分解并重新合項(xiàng)可得:

式中:y0(αeff)為對應(yīng)于非定常分離流的強(qiáng)制函數(shù)。為了建立模型方便，采用如下函數(shù)形式［5］:

2.2 橫向狀態(tài)空間建模原理

建模前首先將飛行器分解成左機(jī)翼、右機(jī)翼及垂直尾翼三個(gè)面。其原因就是為了使用各個(gè)面上的當(dāng)?shù)赜羌敖撬俣确謩e建立相應(yīng)的法向氣動(dòng)力狀態(tài)空間模型，最后得到橫向氣動(dòng)力矩。由理論力學(xué)定義可得某一點(diǎn)的力矩:

由定義可知，滾轉(zhuǎn)力矩是飛行器氣動(dòng)力關(guān)于體軸系下x軸的力矩。對于具體的問題，由于飛行器橫向的對稱性，模型的復(fù)雜性就大大減小。通過這種分解的思想，就可以單獨(dú)描述每個(gè)面上的氣動(dòng)力特性。應(yīng)用這種方法對每一個(gè)面上的流場建模時(shí)，模型能夠反映一個(gè)機(jī)翼已經(jīng)分離而另外一個(gè)機(jī)翼卻很少分離的物理現(xiàn)象。

式中:p為繞飛機(jī)體軸系x軸的角速度;yp，zp分別為機(jī)翼和垂尾上任意一點(diǎn)在體軸系下y方向和z方向的坐標(biāo)值。對于翼面任意一點(diǎn)處當(dāng)?shù)赜请S時(shí)間的變化率，直接對上兩式求導(dǎo)即可。

在確定了升力面分解的思路和角度轉(zhuǎn)換關(guān)系后，下面推導(dǎo)具體的橫向狀態(tài)空間表達(dá)式。對于橫向運(yùn)動(dòng)，各翼面的法向力系數(shù)CN是補(bǔ)充狀態(tài)變量y(t)和輸入變量α(t)及)的泛函。因此，飛機(jī)滾轉(zhuǎn)力矩系數(shù)Cl也是上述變量的泛函。

因?yàn)闈L轉(zhuǎn)力矩系數(shù)是法向力系數(shù)的函數(shù)，所以，僅僅需要對法向力系數(shù)進(jìn)行泰勒展開，而不用直接對滾轉(zhuǎn)力矩系數(shù)進(jìn)行泰勒展開。每個(gè)面上法向力系數(shù)的泰勒展開式如下:

現(xiàn)在再用一次泰勒展開式來推導(dǎo)式(12)中的這些偏導(dǎo)數(shù)與分離狀態(tài)變量y的關(guān)系表達(dá)式。以法向力系數(shù)對迎角的一階偏導(dǎo)數(shù)CNα(y)為例，先關(guān)于常量值進(jìn)行泰勒展開，其中 ^y∈［0，1］:

忽略二階以上的高階項(xiàng)，并重組方程可得:

依此方法，便可得到各翼面處法向力系數(shù)偏導(dǎo)數(shù)的具體表達(dá)式如下:

從式(15)可以發(fā)現(xiàn)參數(shù) ai(i=1，2，3，4)相當(dāng)于常規(guī)穩(wěn)定導(dǎo)數(shù)，而參數(shù)bi和ci則與流動(dòng)分離或渦破裂位置相關(guān)。

滾轉(zhuǎn)力矩系數(shù)的計(jì)算公式為:

式中:ym為機(jī)翼上氣動(dòng)力中心到x軸的距離;zm為垂直尾翼上的氣動(dòng)中心到x軸的距離。它們都由參數(shù)辨識(shí)得到。

至此，就得到了橫向滾轉(zhuǎn)力矩模型。需要說明的是，實(shí)際應(yīng)用中因?yàn)榇刮膊皇且粋€(gè)大的影響因素，并且垂直尾翼附近的流場受機(jī)翼表面流場的影響極大，因此將垂尾法向力系數(shù)對迎角和迎角變化率的偏導(dǎo)數(shù)當(dāng)作待辨識(shí)的常數(shù)處理。

2.3 辨識(shí)準(zhǔn)則及參數(shù)辨識(shí)

數(shù)學(xué)模型確定后，就需要建立辨識(shí)準(zhǔn)則并運(yùn)用滾轉(zhuǎn)振蕩試驗(yàn)數(shù)據(jù)，通過相關(guān)辨識(shí)算法確定模型中的待定參數(shù)。本文采用最小均方差準(zhǔn)則作為辨識(shí)準(zhǔn)則。傳統(tǒng)的參數(shù)辨識(shí)算法如牛頓-拉夫遜法等進(jìn)行求解時(shí)是基于梯度搜索方法進(jìn)行的，在迭代過程中容易陷入局部最小值，現(xiàn)介紹一種具有較好全局搜索功能的粒子群算法［6］。

假設(shè)在一個(gè)D維的搜索空間中，有n個(gè)粒子組成的種群 X=(X1，X2，…，Xn)，其中第 i個(gè)粒子表示為一個(gè) D 維的向量Xi=［xi1，xi2，…，xiD］T，且表示第i個(gè)粒子在D維搜索空間內(nèi)的位置，亦代表問題的一個(gè)潛在解。根據(jù)辨識(shí)準(zhǔn)則計(jì)算出每個(gè)粒子位置Xi對應(yīng)的適應(yīng)度值。第 i個(gè)粒子的速度為Vi=［Vi1，Vi2，…，ViD］T，其個(gè)體極值為 Pi=［Pi1，Pi2，…，PiD］T，種群的全局極值為 Pg=［Pg1，Pg2，…，PgD］T。在每一次迭代過程中，粒子通過個(gè)體極值和全局極值更新自身的速度和位置:

式中:Vid為粒子的速度;ω為慣性權(quán)重;i=1，2，…，n;d=1，2，…，D;k為當(dāng)前迭代次數(shù);c1和 c2為非負(fù)的常數(shù)，稱為加速因子;r1和r2為分布于［0，1］之間的隨機(jī)數(shù)。

3 模型算例及結(jié)果分析

在應(yīng)用大幅滾轉(zhuǎn)振蕩試驗(yàn)數(shù)據(jù)對上述橫向狀態(tài)空間模型進(jìn)行參數(shù)辨識(shí)的過程中，未知參數(shù)有σ，和zm，其中a1表征小迎角時(shí)完全附著流法向力系數(shù)的靜導(dǎo)數(shù)，是已知常數(shù)。應(yīng)用縱向靜態(tài)試驗(yàn)數(shù)據(jù)建立的狀態(tài)空間模型可以辨識(shí)得到參數(shù)σ，α*，b1和c1，其他參數(shù)用動(dòng)態(tài)試驗(yàn)數(shù)據(jù)辨識(shí)得到。需要注意的是，靜態(tài)數(shù)據(jù)辨識(shí)得到的參數(shù)b1，c1將作為動(dòng)態(tài)試驗(yàn)數(shù)據(jù)參數(shù)辨識(shí)的初始輸入，而σ，α*保持不變。

本次建模共使用3組試驗(yàn)數(shù)據(jù)參與辨識(shí)過程，表1給出了樣本的振幅和頻率。

表1 滾轉(zhuǎn)振蕩運(yùn)動(dòng)(θ=20°)Table 1 Roll oscillatory maneuvers(θ=20°)

表中:R1為模型檢驗(yàn)樣本;R2和R3為建模訓(xùn)練樣本。滾轉(zhuǎn)振蕩時(shí)，每周期采集240個(gè)點(diǎn)，建模時(shí)每隔4個(gè)點(diǎn)取1點(diǎn)，總共120個(gè)樣本作為模型輸入。

圖3和圖4給出了模型和訓(xùn)練樣本的滾轉(zhuǎn)力矩系數(shù)曲線，圖5則給出了模型預(yù)測滾轉(zhuǎn)力矩系數(shù)曲線。從模型預(yù)測結(jié)果來看，模型預(yù)測和試驗(yàn)數(shù)據(jù)具有較好的一致性，而且模型能夠體現(xiàn)不同頻率和振幅的影響，表現(xiàn)了橫向振蕩運(yùn)動(dòng)中非定常氣動(dòng)力的遲滯特性及非線性特性。

圖3 模型和訓(xùn)練樣本對比(R2)Fig.3 Comparison between model and test(R2)

圖4 模型和訓(xùn)練樣本對比(R3)Fig.4 Comparison between model and test(R3)

圖5 模型預(yù)測曲線(R1)Fig.5 Model prediction graph(R1)

4 結(jié)束語

本文探索了應(yīng)用升力面分解的思想建立橫向狀態(tài)空間模型的方法，并初步驗(yàn)證了該方法的可行性，為今后研究適用于任意機(jī)動(dòng)的狀態(tài)空間模型提供了新的思路。但在工作中，僅僅用固定俯仰角下不同頻率和振幅的試驗(yàn)數(shù)據(jù)建立橫向狀態(tài)空間模型，而適用于不同俯仰角下滾轉(zhuǎn)振蕩運(yùn)動(dòng)的模型還有待更廣泛的研究和驗(yàn)證。此外，本文建模過程中，把機(jī)翼氣動(dòng)力中心到x軸的距離ym及垂直尾翼上氣動(dòng)中心到x軸的距離zm當(dāng)作待辨識(shí)常數(shù)來處理，以簡化辨識(shí)過程，而實(shí)際模型滾轉(zhuǎn)振蕩時(shí)，這些參數(shù)隨滾轉(zhuǎn)角和滾轉(zhuǎn)角速度而變化，這都是下一步研究工作中需要細(xì)化和解決的難題。

［1］ Tobak M，Schiff L B.On the formulation of the aerodynamic characteristics in aircraft dynamics［R］.NASATRR-456，1976.

［2］ Goman M G，Khrabrov A N.State-space representation of aerodynamic characteristics of an aircraft at high angels of attack［J］.Journal of Aircraft，1994，31(5):1109-1115.

［3］ 史志偉，吳根興.多變量非線性非定常氣動(dòng)力的模糊邏輯模型［J］.空氣動(dòng)力學(xué)學(xué)報(bào)，2001，19(1):103-108.

［4］ Stagg G A.An unsteady aerodynamic model for use in the high angle of attack regime［D］.Virginia:Virginia Polytechnic Institute and State University，1998.

［5］ Fan Yigang.Identification of an unsteady aerodynamic model up to high angle of attack regime［D］.Virginia:Virginia Polytechnic Institute and State University，1997.

［6］ 柯晶，錢積新.應(yīng)用粒子群優(yōu)化的非線性系統(tǒng)辨識(shí)［J］.電路與系統(tǒng)學(xué)報(bào)，2003，8(4):12-16.