含分?jǐn)?shù)階導(dǎo)數(shù)項(xiàng)的隨機(jī)Duffing振子的穩(wěn)態(tài)響應(yīng)分析

孫春艷, 徐 偉

(1．山東大學(xué)(威海)數(shù)學(xué)與統(tǒng)計(jì)學(xué)院，山東 威海 264209； 2.西北工業(yè)大學(xué)應(yīng)用數(shù)學(xué)系，陜西 西安 710072)

含分?jǐn)?shù)階導(dǎo)數(shù)項(xiàng)的隨機(jī)Duffing振子的穩(wěn)態(tài)響應(yīng)分析

孫春艷1, 徐 偉2

(1．山東大學(xué)(威海)數(shù)學(xué)與統(tǒng)計(jì)學(xué)院，山東 威海 264209； 2.西北工業(yè)大學(xué)應(yīng)用數(shù)學(xué)系，陜西 西安 710072)

對一個含有分?jǐn)?shù)階導(dǎo)數(shù)項(xiàng)阻尼的、Gaussian白噪聲激勵下的Duffing振子進(jìn)行了穩(wěn)態(tài)響應(yīng)分析。首先，基于能量平衡理論，運(yùn)用等效線性化方法，計(jì)算等效系統(tǒng)的線性阻尼及自然頻率，建立統(tǒng)計(jì)意義下的等效線性化系統(tǒng)。然后，利用平均法建立隨機(jī)It方程，得到隨機(jī)響應(yīng)的Markovian近似；給出描述振子振幅概率密度函數(shù)演化的Fokker-Planck方程，并得到它的穩(wěn)態(tài)解。進(jìn)一步，對于含有響應(yīng)振幅的等效線性系統(tǒng)，借助由Laplace變換得到的轉(zhuǎn)換函數(shù)，得到原系統(tǒng)的條件功率譜密度，結(jié)合振幅的穩(wěn)態(tài)概率密度作為權(quán)重函數(shù)，給出原系統(tǒng)功率譜密度的估計(jì)，以及響應(yīng)的統(tǒng)計(jì)量的估計(jì)。數(shù)值模擬的結(jié)果說明所提出的功率譜密度的近似解析表達(dá)式是可靠的，它甚至適用于Duffing振子具有強(qiáng)非線性回復(fù)力的情形，因?yàn)樗梢暂^好地表現(xiàn)出功率譜密度共振頻譜加寬及多峰現(xiàn)象的出現(xiàn)。

分?jǐn)?shù)階導(dǎo)數(shù)； 等效線性化法； 隨機(jī)平均法； 條件功率譜密度； 響應(yīng)的功率譜密度估計(jì)

引 言

隨機(jī)響應(yīng)分析是隨機(jī)動力學(xué)研究的重要方面，已有了一些經(jīng)典的方法和結(jié)果[1-3]。對于響應(yīng)的功率譜密度，文獻(xiàn)[4]利用經(jīng)典的等效線性化方法表現(xiàn)出缺陷，它給出的估計(jì)得到了正確的共振頻率，卻高估了共振處的峰值而且低估了頻譜的寬度。一種改進(jìn)的等效線性化方法在文獻(xiàn)[5]中被提出，它與傳統(tǒng)方法的主要區(qū)別在于將等效線性系統(tǒng)的待定參數(shù)設(shè)定為響應(yīng)振幅的函數(shù)，從而建立了一個求解等效的線性剛度的計(jì)算方法。在文獻(xiàn)[6]中，首次出現(xiàn)了術(shù)語“條件功率譜密度”，并通過與隨機(jī)平均法的結(jié)合，以概率密度函數(shù)為權(quán)重函數(shù)，將響應(yīng)的功率譜密度作為條件功率譜密度的加權(quán)和，給出了隨機(jī)響應(yīng)功率譜密度的有效估計(jì)。在文獻(xiàn)[4]中被正式提出后，條件功率譜密度的概念開始被廣泛地應(yīng)用于響應(yīng)功率譜密度的估計(jì)，給出了理想的近似解析結(jié)果，但它的數(shù)學(xué)嚴(yán)密性及合理性的說明卻一直是一個空白，直到在文獻(xiàn)[7]中，Spanos等才給出了嚴(yán)格的證明。

近年來，分?jǐn)?shù)階微積分在結(jié)構(gòu)動力學(xué)及工程力學(xué)領(lǐng)域，越來越多的研究表明分?jǐn)?shù)階導(dǎo)數(shù)可以用較少的參數(shù)來模擬某些黏彈性材料的本構(gòu)關(guān)系[8-10]，在確定性的響應(yīng)分析中，基于特征向量展開，Laplace變換[8-9]和Fourier變換[11]被用來得到精確的解析解。文獻(xiàn)[12]中，平均法被用來分析響應(yīng)的振幅。在隨機(jī)響應(yīng)的情形，已有的方法主要適用于求解弱非線性系統(tǒng)。文獻(xiàn)[13]中，基于廣義諧和函數(shù)的隨機(jī)平均法被用于對高斯白噪聲激勵下的含有分?jǐn)?shù)階導(dǎo)數(shù)的強(qiáng)非線性系統(tǒng)進(jìn)行隨機(jī)響應(yīng)分析。此外，數(shù)值方法方面，已有的結(jié)果如有限差分方法、有限元方法及無網(wǎng)格方法等[14-16]，都一定程度上受到分?jǐn)?shù)階導(dǎo)數(shù)全局依賴性這一本質(zhì)特征對計(jì)算速度的影響。文獻(xiàn)[17]給出的算法在一定程度上減低了分?jǐn)?shù)階導(dǎo)數(shù)項(xiàng)對歷史數(shù)據(jù)的依賴和記憶性，提高了數(shù)值計(jì)算的速度。

本文旨在建立一個對含有分?jǐn)?shù)階導(dǎo)數(shù)項(xiàng)的隨機(jī)Duffing振子進(jìn)行隨機(jī)響應(yīng)分析的方法。利用改進(jìn)的等效線性化方法，首先得到一個含分?jǐn)?shù)階導(dǎo)數(shù)項(xiàng)的線性系統(tǒng)，其中的線性阻尼及自然頻率都是振幅的函數(shù)。將分?jǐn)?shù)階導(dǎo)數(shù)項(xiàng)作為一個弱阻尼項(xiàng)，對得到的系統(tǒng)運(yùn)用平均法，建立隨機(jī)It方程，得到響應(yīng)的Fokker-Planck方程并求得其穩(wěn)態(tài)解。然后，借助Laplace變換，得到等效線性系統(tǒng)的轉(zhuǎn)換函數(shù)及條件功率譜密度，并綜合之前的結(jié)果，對條件功率譜密度利用振幅的概率密度函數(shù)進(jìn)行加權(quán)求和，得到原系統(tǒng)響應(yīng)功率譜密度的估計(jì)。最后，通過數(shù)值模擬的結(jié)果來驗(yàn)證響應(yīng)的穩(wěn)態(tài)概率密度及功率譜密度軌跡的合理性。

1 含分?jǐn)?shù)階導(dǎo)數(shù)項(xiàng)的隨機(jī)Duffing振子

考慮一個高斯白噪聲激勵下的含分?jǐn)?shù)階導(dǎo)數(shù)項(xiàng)的Duffing振子

(1)

(2)

(3)

(4)

并引入?yún)?shù)D=πS0來表征噪聲的強(qiáng)度。

2 等效線性化系統(tǒng)

有別于傳統(tǒng)的等效線性化方法，從能量平衡的角度，改進(jìn)的線性化方法認(rèn)為原振動結(jié)構(gòu)的能耗與等價(jià)系統(tǒng)的能耗相同，以此得到等價(jià)線性系統(tǒng)的待定阻尼系數(shù)，而待定的自然頻率可以通過完全獨(dú)立的其他方法來得到，普遍的方法是將它認(rèn)為是原系統(tǒng)的共振頻率[19]。

將待定的阻尼及剛度系數(shù)作為系統(tǒng)振幅響應(yīng)的函數(shù)，設(shè)等價(jià)線性振動結(jié)構(gòu)為

(5)

(6)

利用廣義諧波平衡技術(shù)，可以得到如下的待定系數(shù)，使誤差(6)達(dá)到均方意義的最小

(7)

至此，可以得到等效的線性系統(tǒng)為

(8)

3 平均法、隨機(jī)響應(yīng)的穩(wěn)態(tài)概率密度

在弱阻尼下隨機(jī)系統(tǒng)的響應(yīng)分析中，平均法可以建立原系統(tǒng)振幅響應(yīng)的Markovian近似，并通過求解描述振幅概率密度函數(shù)演化的Fokker-Planck方程，得到隨機(jī)響應(yīng)的概率分布及統(tǒng)計(jì)特性。對于線性系統(tǒng)(8)，當(dāng)阻尼系數(shù)ξ是一個小量時(shí)，它在原點(diǎn)附近的運(yùn)動可以被看作是周期的，據(jù)此引入變換

(9)

(10)

(11)

這里m(A)和σ(A)分別為待定的漂移和擴(kuò)散系數(shù)，B(t)是單位Wiener過程。

(12)

為計(jì)算m(A)，需要分別計(jì)算兩個確定性平均，對0<α<1，其中第一項(xiàng)有

(13)

其中

(14a)

(14b)

由于A(t)和Θ(t)都是慢變變量，式(14)中

(15)

即有

(16)

對于式(16)中的Fresnel積分，引入漸近積分

(17a)

(17b)

將式(16)及(17)代入L1和L2，得到

(18)

及

(19)

綜上，當(dāng)0<α<1時(shí)，平均后的結(jié)果為

(20)

對于1<α<2的情形，類似的推導(dǎo)得到

(21)

將式(20)及(21)代入計(jì)算公式(12)，結(jié)合G1和G2的表達(dá)式，得到漂移系數(shù)

m(A)=

(22)

及擴(kuò)散系數(shù)

(23)

對于It隨機(jī)微分方程(11)，描述隨機(jī)響應(yīng)的概率密度函數(shù)變化的Fokker-Planck方程為

(24)

在穩(wěn)態(tài)情形，方程(24)的解為

(25)

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

圖1～4給出了振幅響應(yīng)的穩(wěn)態(tài)概率密度，通過不同階數(shù)α及不同非線性程度ε下近似解析結(jié)果與數(shù)值模擬結(jié)果的對比，說明式(25)給出的結(jié)果的合理性。其中數(shù)值模擬采用的參數(shù)分別為：D=0.1，ξ=0.1，ω0=2，其他參數(shù)取值不同，分別在圖中給予說明。從圖1和2中可以看出，當(dāng)分?jǐn)?shù)階階數(shù)α取值為0.5時(shí)，對于強(qiáng)非線性系統(tǒng)情形，即ε=1及ε=3，利用隨機(jī)平均法所得到的穩(wěn)態(tài)概率密度函數(shù)與響應(yīng)的數(shù)值模擬結(jié)果十分貼近。這說明，在0<α<1時(shí)，所得到的穩(wěn)態(tài)概率密度的近似解析結(jié)果是有效的。類似地，圖3和4給出了分?jǐn)?shù)階階數(shù)α取值為1.5時(shí)穩(wěn)態(tài)概率密度解析解與模擬解的對比情況，對ε=1及ε=3兩種情形，近似解析解與數(shù)值模擬解很好的吻合，說明所得到的解析解在1<α<2時(shí)也是有效的。

圖1 α=0.5，ε=1時(shí)振幅的穩(wěn)態(tài)概率密度Fig.1 The stationary probability density function for α=0.5, ε=1

圖2 α=0.5，ε=3時(shí)振幅的穩(wěn)態(tài)概率密度Fig.2 The stationary probability density function for α=0.5, ε=3

圖3 α=1.5，ε=1時(shí)振幅的穩(wěn)態(tài)概率密度Fig.3 The stationary probability density function for α=1.5，ε=1

圖4 α=1.5，ε=3時(shí)振幅的穩(wěn)態(tài)概率密度Fig.4 The stationary probability density function for α=1.5，ε=3

4 功率譜密度估計(jì)

對于得到的等效線性系統(tǒng)(8)，在系統(tǒng)達(dá)到穩(wěn)態(tài)時(shí)，任意給定一個振幅的觀測值A(chǔ)=a，即可得到一個等效系統(tǒng)

(26)

與系統(tǒng)(26)對應(yīng)的轉(zhuǎn)換函數(shù)為

(27)

由此可以得到系統(tǒng)(26)響應(yīng)的功率譜密度

(28)

(29)

隨即可以得到響應(yīng)的統(tǒng)計(jì)特征，如原系統(tǒng)穩(wěn)態(tài)位移響應(yīng)的方差為

(30)

5 數(shù)值結(jié)果

為了驗(yàn)證式(29)給出的響應(yīng)的功率譜密度估計(jì)，將近似的解析結(jié)果與數(shù)值模擬的結(jié)果進(jìn)行對比。圖5～8中，參數(shù)分別為：D=0.1，ξ=0.01，ω0=2，Welch法被運(yùn)用于響應(yīng)的數(shù)值模擬結(jié)果，得到響應(yīng)功率譜密度的Welch法估計(jì)。觀察圖5及6，在分?jǐn)?shù)階階數(shù)α取值為0.5時(shí)，對于ε=1及ε=3，對比式(29)給出的近似解析結(jié)果與數(shù)值結(jié)果，響應(yīng)的功率譜密度的估計(jì)不僅給出了正確的共振位置，也給出了合理的共振頻譜帶寬。說明在0<α<1時(shí)，利用條件功率譜密度所建立的式(29)給出的響應(yīng)的功率譜密度估計(jì)是有效的。同樣地，對比圖7及8，在分?jǐn)?shù)階階數(shù)α取值為1.5時(shí)，對于ε=1及ε=3的情形，對比式(29)給出的近似解析結(jié)果與數(shù)值結(jié)果，響應(yīng)的功率譜密度的估計(jì)都給出了正確的共振位置及合理的共振頻譜帶寬，并有效地表現(xiàn)共振出頻譜加寬及多峰現(xiàn)象的出現(xiàn)。說明在1<α<2時(shí)，利用條件功率譜密度所建立的式(29)給出的響應(yīng)的功率譜密度估計(jì)是有效的。為進(jìn)一步說明響應(yīng)功率譜密度估計(jì)的有效性，利用式(30)所得到的位移響應(yīng)方差的近似解析結(jié)果與模擬結(jié)果進(jìn)行對比，得到圖9，說明所得到的響應(yīng)的功率譜密度估計(jì)可以很好地用于計(jì)算響應(yīng)的統(tǒng)計(jì)特性。

圖5 α=0.5，ε=1時(shí)響應(yīng)的功率譜密度估計(jì)Fig.5 Response power spectral density estimation for α=0.5, ε=1

圖6 α=0.5，ε=3時(shí)響應(yīng)的功率譜密度估計(jì)Fig.6 Response power spectral density estimation for α=0.5, ε=3

圖7 α=1.5，ε=1時(shí)響應(yīng)的功率譜密度估計(jì)Fig.7 Response power spectral density estimation for α=1.5, ε=1

圖8 α=1.5，ε=3時(shí)響應(yīng)的功率譜密度估計(jì)Fig.8 Response power spectral density estimation for α=1.5, ε=3

圖9 位移響應(yīng)的方差Fig.9 Variance for the displacement response

6 結(jié) 論

對于含分?jǐn)?shù)階導(dǎo)數(shù)項(xiàng)的強(qiáng)非線性隨機(jī)Duffing振子，本文給出了一個響應(yīng)功率譜密度估計(jì)的程序。在假設(shè)分?jǐn)?shù)階導(dǎo)數(shù)項(xiàng)表現(xiàn)為結(jié)構(gòu)阻尼的前提下，本文給出的功率譜密度估計(jì)程序是可靠的，即使是在較強(qiáng)的非線性系統(tǒng)情形。首先，利用改進(jìn)的等效線性化法，得到一個與原系統(tǒng)等效的含分?jǐn)?shù)階導(dǎo)數(shù)項(xiàng)的線性隨機(jī)Duffing振子，通過計(jì)算等價(jià)系統(tǒng)的線性阻尼和自然頻率，原系統(tǒng)的非線性回復(fù)力被一個線性的回復(fù)力代替，而且它是振幅響應(yīng)的函數(shù)。然后，隨機(jī)平均法被作用于得到的線性振子，通過系統(tǒng)振幅響應(yīng)的Markovian近似，建立了振幅響應(yīng)的隨機(jī)It方程；寫出描述振幅概率密度函數(shù)演化的Fokker-Planck方程，并給出其穩(wěn)態(tài)解。最后，利用等效線性系統(tǒng)的轉(zhuǎn)換函數(shù)，得到隨機(jī)響應(yīng)的條件功率譜密度，并結(jié)合得到的振幅穩(wěn)態(tài)概率密度函數(shù)，得到響應(yīng)功率譜密度的估計(jì)。數(shù)值結(jié)果表明，文中所給出的近似解析結(jié)果，無論是振幅的穩(wěn)態(tài)概率密度還是響應(yīng)的功率譜密度，對于包括強(qiáng)非線性回復(fù)力情形在內(nèi)的穩(wěn)態(tài)響應(yīng)分析都是有效的，尤其對于強(qiáng)非線性系統(tǒng)，所給出的估計(jì)可以表現(xiàn)出響應(yīng)功率譜密度頻譜的加寬及高次諧波的出現(xiàn)。這在一定程度上說明“改進(jìn)的”等效線性化方法是可靠的，并且，將分?jǐn)?shù)階導(dǎo)數(shù)項(xiàng)作為一個弱阻尼項(xiàng)的假設(shè)是合理的，且基于此建立的隨機(jī)平均法是適用的；條件功率譜密度在含分?jǐn)?shù)階導(dǎo)數(shù)項(xiàng)的系統(tǒng)中仍然適用，說明了它是一個對隨機(jī)響應(yīng)進(jìn)行功率譜密度估計(jì)的有效工具。

[1] Robert J B, Spanos P D. Stochastic averaging: an approximate method of solving random vibration problems[J]. International Journal of Non-linear Mechanics, 1986, 21(2): 111—134.

[2] Robert J B, Spanos P D. Random Vibration and Statistical Linearization[M]. New York: Dover Publications, 2003.

[3] Zhu W Q. Recent developments and applications of the stochastic averaging method in random vibration[J]. Applied Mechanics Reviews, 1996, 49(10S): S72.

[4] Bouc R. The power spectral density of response for a strongly non-linear random oscillator[J]. Journal of Sound and Vibration, 1994, 175(3): 317—331.

[5] Miles R N. An approximate solution for the spectral response of Duffing′s oscillator with random input[J]. Journal of Sound and Vibration, 1989, 132(1): 43—49.

[6] Miles R N. Spectral response of a bilinear oscillator[J]. Journal of Sound and Vibration， 1993, 163(2): 319—326.

[7] Spanos P D, Kougioumtzoglou I A, Soize C. On the determination of the power spectrum of randomly excited oscillators via stochastic averaging: An alternative perspective[J]. Probabilistic Engineering Mechanics, 2011, 26(1): 10—15.

[8] Bagley R L, Torvik P J. Fractional calculus-a different approach to the analysis of viscoelastically damped structures[J]. AIAA Journal, 1983, 21(5): 741—748.

[9] Bagley R L, Torvik P J. Fractional calculus in the transient analysis of viscoelastically damped structures[J]. AIAA Journal, 1985, 23(6): 918—925.

[10]Koeller R C. Application of fractional calculus to the theory of viscoelasticity[J]. ASME Journal of Applied Mechanics, 1984, 51(2): 299—307.

[11]Gaul L, Klein P, Kemple S. Impulse response function of an oscillator with fractional derivative in damping description[J]. Mechanics Research Communications, 1989, 16(5): 297—305.

[12]Wahi P, Chatterjee A. Averaging for oscillations with light fractional order damping[A]. Proceedings of ASME 2003 Design Engineering Technical Conferences and Computers and Information in Engineering Conference[C]． Chicago, 2003: 721—727.

[13]Huang Z L, Jin X L. Response and stability of a SDOF strongly nonlinear stochastic system with light damping modeled by a fractional derivative[J]. Journal of Sound and Vibration, 2009, 319(3-5): 1 121—1 135.

[14]Diethelm K, Ford N J, Freed A D et al. Algorithms for the fractional calculus: A selection of numerical methods[J]. Computer Methods in Applied Mechanics and Engineering, 2005, 194(6-8): 743—773.

[15]Diethelm K, Ford N J, Freed A D. A predictor-corrector approach for the numerical solution of fractional differential equations[J]. Nonlinear Dynamics, 2002, 29(1-4): 3—22.

[16]Yuan L, Agrawal O P. A numerical scheme for dynamic system containing fractional derivatives[J]. Journal of Vibration and Acoustics, 2002, 124(2): 321—324.

[17]Spanos P D, Evangelatos G I. Response of a non-linear system with restoring forces governed by fractional derivatives-time domain simulation and statistical linearization solution[J]. Soil Dynamics and Earthquake Engineering, 2010, 30(9): 811—821.

[18]Podlubny I. Fractional Differential Equations[M]. London: Academic Press, 1999.

[19]Goto H, Iemura H. Linearization techniques for earthquake response of simple hysteretic structures[J]. Proceedings of the Japaneese Society of Civil Engineering, 1973, 212:109—119.

Stationary response analysis for a stochastic Duffing oscillator comprising fractional derivative element

SUNChun-yan1,XUWei2

(1．School of Mathematics and Statistics, Shandong University, Weihai 264209, China;2．Department of Applied Mathematics, Northwestern Polytechnical University, Xi′an 710072, China)

Stationary response is investigated for a Duffing oscillator comprising fractional derivative elements excited by Gaussian white noise in the present paper. Firstly, harmonic balance technique is adopted to form a statistically equivalent linear system. Then, stochastic averaging is applied to the system to obtain a Markovian approximation of the response amplitude, and the associated Fokker-Planck equation and its stationary solution are derived. Furthermore, in virtue of Laplace transform, the transfer function of the equivalent linear system with amplitude-dependent coefficients is derived and it gives the conditional power spectral density, after weighted by the stationary probability density function, estimations of the power spectral density for the response and related statistics are derived. Numerical simulations verify the reliability of the proposed procedure, even for strongly nonlinear oscillators with properties like spectrum broadening and multimodal pattern.

fractional derivative; equivalent linearization; stochastic averaging; conditional power spectral density; response power spectral density estimation

2014-03-26；

2014-09-01

國家自然科學(xué)基金資助項(xiàng)目(11772233)

O322； O324

A

1004-4523(2015)03-0374-07

10.16385/j.cnki.issn.1004-4523.2015.03.006

孫春艷(1984—)，女，博士。電話：13679122401；E-mail: sunchunyan@mail.nwpu.edu.cn