焦萌倩,彭如月,黃文韜,蔣貴榮*

(1.桂林電子科技大學(xué) 數(shù)學(xué)與計算科學(xué)學(xué)院，廣西 桂林 541004；2.廣西師范大學(xué) 數(shù)學(xué)與統(tǒng)計學(xué)院，廣西 桂林 541006)

隨著科學(xué)技術(shù)的飛速發(fā)展，飛行器需要完成各種各樣的高難度動作，機(jī)翼搖滾是常見的橫向不穩(wěn)定現(xiàn)象，因而研究飛行器滾轉(zhuǎn)運(yùn)動已經(jīng)成為當(dāng)下科技發(fā)展的需要。飛行器在實(shí)際飛行中遇到外部環(huán)境以及內(nèi)在動因所帶來的各種隨機(jī)噪聲是不可忽視的。Konstadinopoulos等[1]利用非定常渦格法分析了低速條件下細(xì)長三角翼搖晃的現(xiàn)象，這種方法便于計算，但其前提假定導(dǎo)致應(yīng)用受到很大的限制和約束；Elzebda等[2]比較了細(xì)長三角翼的3種分析模型，發(fā)現(xiàn)僅包含二次項(xiàng)的原始模型不能預(yù)測滾轉(zhuǎn)發(fā)散，并在修正模型基礎(chǔ)上得到運(yùn)動方程解的漸近逼近；劉偉等[3]對靜穩(wěn)定下的細(xì)長三角翼運(yùn)動出現(xiàn)Hopf分叉的臨界條件進(jìn)行了分析，并利用耦合三維非定常Navier-Stokes方程與Euler剛體運(yùn)動方程進(jìn)行數(shù)值模證實(shí)了其理論分析。

隨機(jī)動力學(xué)一直都是研究熱點(diǎn)[4]。Zhu等[5]利用廣義調(diào)和函數(shù)提出了針對單自由度強(qiáng)非線性振子的隨機(jī)平均法，并對Duffing-vanderPol振蕩器在寬帶隨機(jī)參激和外激共同作用下的響應(yīng)進(jìn)行了分析,還利用該方法分析了帶有2種參激的Duffing-vanderPol振蕩器系統(tǒng)的平穩(wěn)概率密度和Hopf分岔；黃志龍等[6]利用廣義諧和函數(shù)的隨機(jī)平均法及Galerkin法研究了受高斯白噪聲激勵帶有時滯反饋控制的單自由度強(qiáng)非線性系統(tǒng)，通過將其轉(zhuǎn)化為不具時滯的系統(tǒng)再利用隨機(jī)平均法求得響應(yīng)的近似瞬態(tài)概率密度；Miles[7]提出了一種近似方法，用于估算Duffing振子在用高斯白噪聲驅(qū)動時的響應(yīng)功率譜密度并發(fā)現(xiàn)隨著隨機(jī)激發(fā)的幅度增加，響應(yīng)譜中的諧振響應(yīng)峰值趨于變寬；黃開嬌等[8]研究了含Lévy噪聲的隨機(jī)捕食-被捕食系統(tǒng)，利用李雅普諾夫函數(shù)等得到系統(tǒng)存在唯一全局正解；仲淑娟等[9]利用李雅普諾夫函數(shù)和多尺度法，針對諧和與隨機(jī)噪聲下的非線性系統(tǒng)的響應(yīng)進(jìn)行了研究。

近幾年，很多學(xué)者也開始注重隨機(jī)理論在實(shí)際應(yīng)用中的研究[10-12]。譚建國等[13]討論了模型在隨機(jī)參數(shù)激勵下的穩(wěn)定性與可靠性，在原有二元機(jī)翼模型的基礎(chǔ)上加入隨機(jī)因素，研究了噪聲強(qiáng)度對系統(tǒng)的可靠性的影響，但其僅考慮了隨機(jī)氣流速度的影響，涉及到的隨機(jī)因素還不夠全面；葛根等[14]先運(yùn)用Galerkin變分法將矩形薄板的受面內(nèi)隨機(jī)激勵的振動模型簡化為常微分非線性動力學(xué)方程，并利用擬不可積Hamilton隨機(jī)平均及伊藤方程微分法則等理論將方程簡化為一維隨機(jī)微分方程進(jìn)行研究，通過穩(wěn)態(tài)概率密度函數(shù)的形狀分析了系統(tǒng)的隨機(jī)分岔現(xiàn)象；覃英華等[15]采用數(shù)值模擬的方法分析了噪聲對電力網(wǎng)絡(luò)穩(wěn)定性的影響，隨機(jī)噪聲使得電力網(wǎng)絡(luò)最終達(dá)到崩潰，因此在實(shí)際應(yīng)用中噪聲不容忽視。

隨機(jī)系統(tǒng)的平穩(wěn)概率密度的計算難度很大，一般情況下很難得到精確表達(dá)式。在已有文獻(xiàn)中，大部分學(xué)者在研究三角翼滾轉(zhuǎn)運(yùn)動時考慮參激或者外激，為了理論分析系統(tǒng)的隨機(jī)響應(yīng)考慮特殊的微分方程。本文在更為復(fù)雜的微分方程的基礎(chǔ)上，同時引入外激和速度參激，利用耗散能量平衡法、幅值包線隨機(jī)平均法、能量包線隨機(jī)平均法求得三角翼滾轉(zhuǎn)運(yùn)動系統(tǒng)近似平穩(wěn)概率密度的解析解，并進(jìn)行數(shù)值模擬。

本文內(nèi)容安排如下：第1章建立隨機(jī)三角翼滾轉(zhuǎn)運(yùn)動模型；第2章運(yùn)用耗散能量平衡法，通過求系統(tǒng)的概率勢直接得到近似概率密度；第3章通過進(jìn)行vanderPol變換，運(yùn)用幅值包線隨機(jī)平均法求得系統(tǒng)近似概率密度；第4章運(yùn)用能量包線隨機(jī)平均法求得近似解；第5章給出數(shù)值模擬和結(jié)論。

1 系統(tǒng)模型

考慮單自由度滾轉(zhuǎn)運(yùn)動方程[16]

(1)

(2)

將式(2)代入式(1)可得

(3)

式中ω2=D-ca1，α=-ca2，β=-ca3，γ=-ca4，η=-ca5。

在模型(3)的基礎(chǔ)上考慮隨機(jī)因素，加入噪聲分別為速度乘性噪聲和加性噪聲。當(dāng)噪聲均為高斯白噪聲時，模型為

(4)

式中W1(t)和W2(t)是相互獨(dú)立高斯白噪聲，其功率譜密度與自相關(guān)函數(shù)分別為

ΦWiWi(ω0)=Kii,RWiWi(τ)=2πKiiδ(τ),i=1,2,

式中：δ(τ)為狄拉克函數(shù)；Kii為高斯白噪聲的強(qiáng)度。

2 耗散能量平衡法

非線性系統(tǒng)只有在特定條件下才可以求得其精確解，但實(shí)際很難滿足限定條件，因此，人們發(fā)展了一系列求解非線性系統(tǒng)近似解的方法。耗散能量平衡法是等效非線性系統(tǒng)法的一種，其可以直接找出系統(tǒng)的概率勢，減少誤差的產(chǎn)生。為應(yīng)用耗散能量平衡法，引入狀態(tài)變量

將式(4)寫成單自由度非線性隨機(jī)系統(tǒng)標(biāo)準(zhǔn)形式

(5)

按Wong-Zakai修正項(xiàng)識別附加恢復(fù)力u*(X1)與附加阻尼力h*(X1,X2)的關(guān)系如下

由此可得

u*(X1)=0,h*(X1,X2)=πK22X2。

(6)

系統(tǒng)(4)恢復(fù)力是線性的，可做如下變換

(7)

由式(5)、(6)及變換(7)可得其概率勢函數(shù)導(dǎo)數(shù)為

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

p(x1,x2)=Cexp[-ψ(λ)]=

(8)

3 幅值包線隨機(jī)平均法

幅值能量平衡法主要應(yīng)用于恢復(fù)力是線性的系統(tǒng)，該方法主要任務(wù)是進(jìn)行vanderPol變換，將系統(tǒng)變?yōu)槭苄》蔷€性阻尼與弱激勵擾動的擬線性振子的幅值與相位，并對其在擬周期上進(jìn)行時間平均。對系統(tǒng)(4)進(jìn)行如下vanderPol變換[17]

(9)

(10)

(11)

將A(t)看作X1(t)，φ(t)看作X2(t)，利用式(10)、(11)可得下列標(biāo)準(zhǔn)形式：

(12)

(13)

(14)

(15)

(16)

A(t)和φ(t)構(gòu)成一個馬爾可夫擴(kuò)散過程，且系統(tǒng)以2π/ω為擬周期運(yùn)動，故可對式(13)～(16)在一個擬周期上做如下時間平均

并運(yùn)用伊藤微分規(guī)則得幅值的伊藤方程。平均后幅值A(chǔ)(t)的方程不含相位φ(t)，故得到幅值的光滑伊藤方程為

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

(17)

式中：

(18)

(19)

對應(yīng)的FPK方程為

為求解該FPK方程，需先求解A(t)的與時間無關(guān)的平穩(wěn)概率密度p(A)。易知該平穩(wěn)概率密度滿足下列簡化的FPK方程

(20)

對式(20)積分得

(21)

由式(21)可知概率G(A)在任何處均為常數(shù)。為確保平穩(wěn)概率密度p(A)存在，對有限或無窮邊界最常見的情形是概率流為零，故可從式(21)得簡化的解

(22)

(23)

4 能量包線隨機(jī)平均法

能量包線隨機(jī)平均也稱為擬保守平均，首先考慮系統(tǒng)的無阻尼自由振動方程為

(24)

勢能與總能量分別為

(25)

周期為

(26)

(27)

由于激勵是高斯白噪聲,可得系統(tǒng)漂移系數(shù)和擴(kuò)散系數(shù)

(28)

(29)

時間平均按式(30)進(jìn)行

(30)

于是，支配能量Λ(t)的伊藤方程為

由此可得Λ(t)的平穩(wěn)概率

(31)

式中C2為積分常數(shù)，且

5 數(shù)值模擬和結(jié)論

在系統(tǒng)(4)中，取ω=0.4,α=0.6,β=1.2,γ=0.8,η=0.49,K11=0.4,K22=0.3。系統(tǒng)(4)過初始點(diǎn)(1,2)的一個樣本解的相圖和時間序列如圖1所示。

圖1 系統(tǒng)(4)在ω=0.4,α=0.6,β=1.2,γ=0.8,η=0.49,K11=0.4,K22=0.3下的一個樣本解的相圖和時間序列Fig.1 Phase diagram and time sequence of a sample solution of system (4)with ω=0.4,α=0.6,β=1.2,γ=0.8,η=0.49,K11=0.4,K22=0.3

利用耗散能量平衡法、幅值包線隨機(jī)平均法、能量包線隨機(jī)平均法求得系統(tǒng)(4)的近似平穩(wěn)概率密度表達(dá)式分別為式(8)、(23)和(31)，分別如圖2(a)、(b)、(c)所示。

圖2 系統(tǒng)(4)在ω=0.4,α=0.6,β=1.2,γ=0.8,η=0.49,K11=0.4,K22=0.3下的平穩(wěn)概率密度Fig.2 Approximate probability density of system (4) withω=0.4,α=0.6,β=1.2,γ=0.8,η=0.49,K11=0.4,K22=0.3

為比較每個方法得到的平穩(wěn)概率密度的精確性，現(xiàn)在計算E(φ2)。根據(jù)近似平穩(wěn)概率密度(8)、(23) 和(31)計算的E(φ2)的圖象在圖3中分別用M-1、M-2、M-3表示。M-4是對系統(tǒng)(4)的φ2的蒙特卡羅模擬結(jié)果，其中樣本容量為1萬個。從圖3(a)可以看出，耗散能量平衡法得到的平穩(wěn)概率密度誤差最大；圖3(b)可以看出，能量包線隨機(jī)平均法得到的平穩(wěn)概率密度最精確。

圖3 系統(tǒng)(4)在ω=0.4,α=0.6,β=1.2,γ=0.8,η=0.49,K22=0.3和K11∈(0.3,0.6)下的E(φ2)Fig.3 E(φ2) of system (4) with K11∈(0.3,0.6) and ω=0.4,α=0.6,β=1.2,γ=0.8,η=0.49,K22=0.3

本文在更為復(fù)雜的描述三角翼滾轉(zhuǎn)運(yùn)動的微分方程的基礎(chǔ)上，同時引入外激和速度參激，利用耗散能量平衡法、幅值包線隨機(jī)平均法、能量包線隨機(jī)平均法求得了三角翼滾轉(zhuǎn)運(yùn)動系統(tǒng)近似平穩(wěn)概率密度的解析解，通過數(shù)值結(jié)果，能量包線隨機(jī)平均法得到的平穩(wěn)概率密度最精確。所得結(jié)果可以為三角翼飛行器滾轉(zhuǎn)運(yùn)動的監(jiān)控策略的制定等提供理論依據(jù)。