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

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

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

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

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

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

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

1 系統(tǒng)模型

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

(1)

(2)

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

(3)

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

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

(4)

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

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

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

2 耗散能量平衡法

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

將式(4)寫成單自由度非線性隨機系統(tǒng)標準形式

(5)

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

由此可得

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

(6)

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

(7)

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

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

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

(8)

3 幅值包線隨機平均法

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

(9)

(10)

(11)

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

(12)

(13)

(14)

(15)

(16)

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

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

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

(17)

式中：

(18)

(19)

對應的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 能量包線隨機平均法

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

(24)

勢能與總能量分別為

(25)

周期為

(26)

(27)

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

(28)

(29)

時間平均按式(30)進行

(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)過初始點(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

利用耗散能量平衡法、幅值包線隨機平均法、能量包線隨機平均法求得系統(tǒng)(4)的近似平穩(wěn)概率密度表達式分別為式(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)可以看出，能量包線隨機平均法得到的平穩(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

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