張永芳，蔣貴榮

(桂林電子科技大學(xué) 數(shù)學(xué)與計(jì)算科學(xué)學(xué)院，廣西 桂林 541004)

三角翼飛機(jī)(delta-wingairplane)是機(jī)翼前緣后掠、后緣基本平直、半翼俯視平面形狀為三角形的飛行器。三角翼機(jī)型有著廣泛的應(yīng)用，如殲-8、米格-21、蘇-15殲擊機(jī)的平尾式三角翼飛機(jī)，“幻影”Ⅲ型殲擊機(jī)和“協(xié)和”式超音速客機(jī)等的無平尾式三角翼飛機(jī)。對于戰(zhàn)斗機(jī)而言，三角翼可以加強(qiáng)結(jié)構(gòu)和氣動(dòng)穩(wěn)定性，從而提高生存率。

當(dāng)飛機(jī)的攻角超過臨界角造成失速時(shí)，控制系統(tǒng)的累計(jì)誤差和外部環(huán)境等都會(huì)造成飛機(jī)失去平穩(wěn)，而飛機(jī)的橫向運(yùn)動(dòng)會(huì)導(dǎo)致其失去平衡性與穩(wěn)定性。飛機(jī)的穩(wěn)定性是指在飛行中受外力干擾后不需要操作人員的干預(yù)，靠自身特性恢復(fù)原來狀態(tài)的能力。飛行器在飛行的過程中，往往電路系統(tǒng)、自動(dòng)控制系統(tǒng)等內(nèi)部因素會(huì)受到隨機(jī)干擾或風(fēng)力、濕度、氣壓等外部因素影響飛行器的平穩(wěn)性，所以在分析飛行器滾轉(zhuǎn)運(yùn)動(dòng)時(shí)，應(yīng)該考慮這些隨機(jī)因素。

Elzebda等[1]對細(xì)長三角翼的亞音速機(jī)翼滾轉(zhuǎn)現(xiàn)象的3種模型進(jìn)行了數(shù)值模擬，給出了修正后的運(yùn)動(dòng)方程的漸近逼近解，但未討論系統(tǒng)加入隨機(jī)的情況。Kori等[2]提出了一種基于擴(kuò)展?fàn)顟B(tài)觀測器(ESO)技術(shù)方法，對細(xì)長三角翼的滾轉(zhuǎn)運(yùn)動(dòng)進(jìn)行魯棒控制，Chavatzopoulos等[3]利用不完整的神經(jīng)網(wǎng)絡(luò)(PNN)來識別三角翼飛機(jī)的滾轉(zhuǎn)運(yùn)動(dòng)，他們所討論的是未考慮隨機(jī)因素的確定性方程。Luo等[4]利用Beecham-Titchener的平均方法，得到了滾轉(zhuǎn)運(yùn)動(dòng)抑制的最優(yōu)控制輸入。陳永亮[5]提出非線性動(dòng)力學(xué)特性分析與控制律設(shè)計(jì)的一體化設(shè)計(jì)思想，應(yīng)用于飛機(jī)大迎角飛行時(shí)機(jī)翼搖晃抑制，揭示機(jī)翼搖晃的全局運(yùn)動(dòng)特性，但控制律過于復(fù)雜，并且抑制結(jié)果是在理想狀態(tài)下得出的。李春濤等[6]設(shè)計(jì)以滾轉(zhuǎn)角速率為內(nèi)回路的滾轉(zhuǎn)角控制律，應(yīng)用魯棒伺服LQR優(yōu)化方法，得到滾轉(zhuǎn)角控制律參數(shù)，驗(yàn)證以滾轉(zhuǎn)角速率為主控變量的控制器具有抗干擾能力，但其控制律過于復(fù)雜。楊利紅[7]對飛行器縱向、橫側(cè)向通道和速度通道進(jìn)行控制仿真設(shè)計(jì)，得到相對簡單的滾轉(zhuǎn)角控制律，但無實(shí)例分析。徐孝誠等[8]用MSC/NASTRAN軟件對再入飛行器復(fù)雜結(jié)構(gòu)軸向隨機(jī)振動(dòng)試驗(yàn)和橫向隨機(jī)振動(dòng)試驗(yàn)進(jìn)行了響應(yīng)分析，吳衛(wèi)國等[9]對飛機(jī)起落架動(dòng)力學(xué)模型進(jìn)行求解與響應(yīng)分析，但他們都未對三角翼飛行器滾轉(zhuǎn)運(yùn)動(dòng)隨機(jī)響應(yīng)進(jìn)行分析。

為此，考慮隨機(jī)因素的三角翼飛行器滾轉(zhuǎn)運(yùn)動(dòng)單自由度系統(tǒng)的隨機(jī)響應(yīng)與控制問題，分析其穩(wěn)定性。

1 飛行器滾轉(zhuǎn)運(yùn)動(dòng)模型

對于三角翼飛機(jī)單自由度滾轉(zhuǎn)運(yùn)動(dòng)，滾轉(zhuǎn)運(yùn)動(dòng)方程為

(1)

(2)

將式(2)代入系統(tǒng)(1)，得

(3)

其中ω2=-ca1,μ1=ca2-D,b1=ca3,μ2=ca4,b2=ca5,μ3=c6a,ai(i=1，2，…，6)是關(guān)于迎角α的函數(shù)，參數(shù)值由表1給出。

表1 系統(tǒng)(3)中的參數(shù)

2 隨機(jī)響應(yīng)分析

考慮高斯白噪聲下三角翼飛行器滾轉(zhuǎn)運(yùn)動(dòng)，其運(yùn)動(dòng)方程為

(4)

其中W(t)為譜密度為K的高斯白噪聲。

隨機(jī)系統(tǒng)(4)的恢復(fù)力u(φ)=ω2φ-b1φ3,其為非線性恢復(fù)力。由于幅值包線隨機(jī)平均要求系統(tǒng)的恢復(fù)力是線性的，系統(tǒng)(4)不能用幅值包線隨機(jī)平均方法討論隨機(jī)響應(yīng)和穩(wěn)定性。采用極坐標(biāo)方法[10]處理系統(tǒng)(4)，考慮極坐標(biāo)變換：

x1=rcosθ,x2=rsinθ,

得到

b1r3sinθcos3θ+μ2r3sin2θcos2θ+

b2r3sin3θcosθ+μ3r3sin4θ+

rsin2θW(t)=f1(r,θ)+g1W(t),

b1r3cos4θ+μ2r3sinθcos3θ+

b2r3sin2θcos2θ+μ3r3sin3θcosθ-

sin2θ+sinθcosθW(t)=f2(r,θ)+g2W(t),

幅值r(t)與相位θ(t)構(gòu)成一個(gè)矢量馬爾可夫擴(kuò)散過程。應(yīng)用隨機(jī)時(shí)間平均法，上式中的fi與gi(i=1，2)是以2π為周期的時(shí)間t的周期函數(shù)，可在一個(gè)周期上進(jìn)行時(shí)間平均，

高斯白噪聲激勵(lì)下系統(tǒng)的漂移系數(shù)m(r)、擴(kuò)散系數(shù)σ(r)為：

σ2(r)=2πK(〈g1g1〉t+〈g1g2〉t)。

其中：

則

(5)

(6)

由于平均后的幅值r(t)方程不含θ(t)，光滑后的幅值r(t)為馬爾可夫擴(kuò)散過程，受伊藤方程支配，

(7)

其中B(t)為單位維納過程，系統(tǒng)(7)為一維伊藤方程。幅值r(t)的平穩(wěn)概率密度為

其中常數(shù)C為歸一化系數(shù)。

系統(tǒng)(7)對應(yīng)的線性方程為

(8)

令Y(t)=lnr(t)，用伊藤微分規(guī)則得到Y(jié)(t)的伊藤方程，

方程的解為

(9)

根據(jù)李雅普諾夫指數(shù)定義，可得到系統(tǒng)(9)的李雅普諾夫指數(shù)為

當(dāng)λ>0，μ1/2+πK/4>0時(shí)，隨機(jī)微分方程(7)對應(yīng)的線性方程的平凡解局部不穩(wěn)定，因此隨機(jī)系統(tǒng)(4)在平衡點(diǎn)(0，0)處局部不穩(wěn)定。

定理1假設(shè)下列條件成立：

其中μ1>0，譜密度K為非負(fù)，則隨機(jī)系統(tǒng)(4)在平衡點(diǎn)(0，0)處局部不穩(wěn)定。

由表1知，當(dāng)飛行器的迎角α>19.5°，μ1>0,因而λ>0，則隨機(jī)系統(tǒng)(4)在平衡點(diǎn)(0，0)處局部不穩(wěn)定。隨著時(shí)間t的增加，飛行器滾轉(zhuǎn)運(yùn)動(dòng)將無法保持穩(wěn)定。

3 滾轉(zhuǎn)運(yùn)動(dòng)的控制

在無人操作的情況下，飛行器自身調(diào)節(jié)達(dá)到平穩(wěn)狀態(tài)，這就需要在飛行器上加入滾轉(zhuǎn)角的控制器。首先飛行器通過測量裝置測出實(shí)際滾轉(zhuǎn)角度，然后將所測信息傳輸給控制器，接著控制器對接收測量信息進(jìn)行處理，當(dāng)滾轉(zhuǎn)角大于期望角度范圍時(shí)，控制器將控制信號傳輸給副翼舵機(jī)系統(tǒng)，最后接到指令的副翼舵機(jī)系統(tǒng)調(diào)節(jié)副翼狀態(tài)，達(dá)到減小滾轉(zhuǎn)角的目的。采用的控制律為

其中：Kp為滾轉(zhuǎn)角傳感器信號；Kd為滾轉(zhuǎn)角速率傳感器信號；φh期望滾轉(zhuǎn)角。圖1為控制滾轉(zhuǎn)角的流程圖，其中eh為期望信號偏差，eh=φh-φ。

圖1 滾轉(zhuǎn)角控制流程圖

在飛行器橫向運(yùn)動(dòng)不穩(wěn)定的狀態(tài)下，加入控制律得到

(10)

采用改進(jìn)隨機(jī)平均方法[11]，求系統(tǒng)(10)的一維幅值伊藤方程。

系統(tǒng)(10)的無阻尼自由振動(dòng)為

(11)

g(φ)=ω2φ-b1φ3-Kpφ,

(12)

(13)

(14)

(15)

(16)

其中ν(A,t)為系統(tǒng)瞬時(shí)頻率。將方程(16)代入方程(14)可得

其中A為系統(tǒng)位移的幅值。將方程(16)代入方程(15)，得到總能量EH，

幅值A(chǔ)與總能量EH的關(guān)系為

G(A)=G(-A)=EH。

通過能量平衡法獲得系統(tǒng)的平均頻率ω(A)，代替該瞬時(shí)頻率ν(A,t)。當(dāng)θ(t)=0時(shí)，

(17)

將平均頻率ω(A)代替瞬時(shí)頻率ν(A,t)代入方程(16)，然后將所得結(jié)果代入方程(15)，得到相關(guān)Hamilton系數(shù)EH,θ=ωt。因此，可得到殘余項(xiàng)

R(t)=EH,θ=ωt-EH,θ=0,

(18)

考慮方程(16)是系統(tǒng)(10)的一個(gè)近似解，所以R不可能總為0值。

當(dāng)ωt=π/4時(shí)，

(19)

方程(19)與式(17)相減，得到余項(xiàng)式

(20)

考慮輸入總能量與輸出總能量近似相等，令

(21)

將方程(17)、(19)、(20)代入方程(21)進(jìn)行求解，可得到其等效的平均頻率ω(A)，

將其代替瞬時(shí)頻率ν(A,t),得到近似方程解

θ(t)=ω(A)t+φ(t),

其中φ(t)為相位。

通過對系統(tǒng)(10)的完全轉(zhuǎn)換化解，可以得到如下方程：

其中：

μ2A3ω(A)2sin2θcos2θ-

b2A3ω(A)3sin3θcosθ+μ3A3ω(A)4sin4θ+

Kpφrω(A)sinθ+KdAω(A)2sin2θ],

μ2A3ω(A)2sinθcos3θ-

b2A3ω(A)3sin2θcos2θ+

μ3A3ω(A)4sin3θcosθ+Kpφrω(A)cosθ+

KdAω(A)2sinθcosθ],

實(shí)施控制后的幅值伊藤方程為

dA=m(A)dt+σ(A)dB(t),

(22)

其中漂移系數(shù)m(A)與擴(kuò)散系數(shù)σ(A)為:

系統(tǒng)(22)的平穩(wěn)概率密度為

[4(ω2-Kp)-3b1A2]δ1Aδ2

其中：

系統(tǒng)(22)對應(yīng)的線性方程為

(23)

令Y(t)=ln[A(t)]，利用伊藤微分規(guī)則得到Y(jié)(t)的伊藤方程：

則其解為

根據(jù)李雅普諾夫指數(shù)定義，可得到系統(tǒng)(15)的李雅普諾夫指數(shù)

當(dāng)λ<0，Kd>μ1+πK/2時(shí)，隨機(jī)微分方程(7)對應(yīng)的線性方程的平凡解局部穩(wěn)定，因此，隨機(jī)系統(tǒng)(4)在平衡點(diǎn)(0，0)處局部穩(wěn)定。

定理2假設(shè)下列條件成立：

其中μ1為正值，譜密度K為非負(fù)，Kd為遠(yuǎn)大于μ1的正值，則隨機(jī)系統(tǒng)(4)在平衡點(diǎn)(0，0)處局部穩(wěn)定。

當(dāng)α>19.5°時(shí)，有μ1>0，當(dāng)Kd>μ1+πK/2時(shí)，系統(tǒng)(10)的李雅普諾夫指數(shù)

則隨機(jī)系統(tǒng)(10)在平衡點(diǎn)(0，0)處局部穩(wěn)定，飛行器在受到隨機(jī)干擾的情況下，通過控制器的調(diào)節(jié)，飛行器的滾轉(zhuǎn)運(yùn)動(dòng)由原來的不穩(wěn)定狀態(tài)變?yōu)榉€(wěn)定狀態(tài)。

4 結(jié)束語

利用隨機(jī)微分方程，建立具有速度參激的三角翼飛行器滾轉(zhuǎn)運(yùn)動(dòng)模型，采用極坐標(biāo)法、時(shí)間平均法與伊藤微分規(guī)則將隨機(jī)系統(tǒng)化為一維伊藤方程，得到系統(tǒng)的平穩(wěn)概率密度和隨機(jī)系統(tǒng)不穩(wěn)定的充分條件。建立加入控制律的滾轉(zhuǎn)運(yùn)動(dòng)新模型，采用改進(jìn)時(shí)間平均法與伊藤微分規(guī)則，將系統(tǒng)化為一維伊藤方程，得到系統(tǒng)穩(wěn)定的充分條件，為三角翼飛行器滾轉(zhuǎn)運(yùn)動(dòng)的控制提供理論依據(jù)。在此條件下，對于系統(tǒng)可能出現(xiàn)的分岔行為，有待進(jìn)一步研究。