韓群 徐偉

(1.華中農(nóng)業(yè)大學理學院,武漢 430070) (2.西北工業(yè)大學應用數(shù)學系,西安 710072)

含指數(shù)積分型阻尼和周期激勵系統(tǒng)的穩(wěn)態(tài)隨機響應分析*

韓群1?徐偉2

(1.華中農(nóng)業(yè)大學理學院,武漢 430070) (2.西北工業(yè)大學應用數(shù)學系,西安 710072)

運用基于短時高斯逼近的廣義胞映射方法,研究了含指數(shù)積分型非粘性阻尼和周期激勵系統(tǒng)在高斯白噪聲作用下的穩(wěn)態(tài)響應.首先介紹了方法的實施過程,并推導了系統(tǒng)的矩方程.然后給出了系統(tǒng)的穩(wěn)態(tài)概率密度函數(shù),分析了阻尼系數(shù)和松弛參數(shù)對穩(wěn)態(tài)響應的影響,并通過直接Monte Carlo模擬的結果驗證了廣義胞映射方法的有效性.

廣義胞映射, 短時高斯逼近, 指數(shù)積分型阻尼, 周期激勵, 穩(wěn)態(tài)響應

引言

廣義胞映射方法已被廣泛用來研究非線性動力系統(tǒng)的隨機響應分析[1-6].傳統(tǒng)廣義胞映射在計算狀態(tài)胞之間一步轉移概率時,通常借助Monte Carlo數(shù)值模擬,對于含周期激勵的系統(tǒng),在考慮一步轉移概率時,轉移時間長度直接取為系統(tǒng)的周期,雖然這樣在形式上顯得很簡單,但是運算過程卻相當費時.因此,該方法在含周期激勵高維系統(tǒng)中的應用顯得非常困難.Sun和Hsu[7]提出一種短時高斯逼近方法構造FPK方程的短時解,從而計算了高斯白噪聲激勵下典型自治系統(tǒng)的一步轉移概率矩陣,該方法顯著提高廣義胞映射方法的效率,在應用中體現(xiàn)出了較大的優(yōu)勢.最近,Han等[8]將基于短時高斯逼近的廣義胞映射方法推廣,用于研究周期和高斯白噪聲激勵下非線性系統(tǒng)的響應概率密度函數(shù),并研究了一類光滑非連續(xù)(SD)振子的瞬態(tài)和穩(wěn)態(tài)響應行為.

指數(shù)積分型非粘性阻尼是一種復雜的阻尼模型[9,10],與經(jīng)典粘性阻尼不同,它不僅與當前的瞬時速度有關,還依賴于速度的歷史狀態(tài).含指數(shù)積分型非粘性阻尼系統(tǒng)的隨機響應研究中,已有一些成果[11,12],但是這些都沒有考慮周期激勵的作用.在處理含這類阻尼的系統(tǒng)時,通常將阻尼看作是一個變量,由此將系統(tǒng)擴展成三維系統(tǒng),若考慮周期和高斯白噪聲共同激勵下該系統(tǒng)的隨機響應分析,必然會對傳統(tǒng)的廣義胞映射方法的應用提出更大的挑戰(zhàn).本文基于文獻[8]中廣義胞映射方法,討論了含指數(shù)積分型非粘性阻尼和周期激勵系統(tǒng)在高斯白噪聲作用下的穩(wěn)態(tài)響應分析.首先簡單介紹含周期激勵系統(tǒng)隨機響應分析的廣義胞映射方法,以及短時高斯逼近在一個周期上的實施策略.然后針對一類含指數(shù)積分型非粘性阻尼的系統(tǒng),先將其擴成三維系統(tǒng),給出了系統(tǒng)的矩方程.最后計算了系統(tǒng)的穩(wěn)態(tài)響應概率密度函數(shù),并討論了阻尼有關參數(shù)對穩(wěn)態(tài)響應的影響.

1 基于短時高斯逼近的廣義胞映射求解過程

考慮含周期和高斯白噪聲激勵的N維非線性隨機動力系統(tǒng),其對應的It隨機微分方程為：

dX(t)=f(X,t)dt+σ(X)dB(t)

(1)

其中X(t)=[X1(t),X2(t),…,XN(t)]T∈RN(上標T表示轉置)是一個N維Markov隨機過程,代表系統(tǒng)的響應.f(X,t)是一個N×1的漂移向量函數(shù),其中包含周期激勵,即滿足f(X,t)=f(X,t+T),T是周期.σ(X)是N×L維耗散矩陣.B(t)=[B1(t),B2(t),…,BL(t)]T∈RL是一個標準的L維向量維納過程,滿足：

E[dB(t)]=0,

(2)

其中I是一個L×L的單位矩陣.維納過程B(t)的形式導數(shù)即為L維的標準高斯白噪聲.

1.1 廣義胞映射方法

n=0,1,2,…

(3)

在廣義胞映射方法的實現(xiàn)過程中,連續(xù)的狀態(tài)空間RN被轉化成離散的胞狀態(tài)空間ZN.考慮一個感興趣的有界區(qū)域Ω,并將之劃分為Nc個常規(guī)胞,它們的尺寸大小相同,假設為h1×h2×…×hN,其中hi是每個胞在xi方向上的長度.所有的常規(guī)胞依次編號為1到Nc, 區(qū)域Ω以外的所有部分被看作是一個陷胞,編號為0. 如果pi(n)表示t=nT時系統(tǒng)的響應位于第i號胞的概率,而qji表示以pi(0)=1為初始條件,系統(tǒng)的響應在T時刻位于第j號胞的轉移概率,那么有：

pi(n)=∫Cip(x,nT)dx

(4)

對于一個常規(guī)胞i,qji可由如下公式計算得到：

(5)

(6)

該方程確定了一個有限Markov鏈的演化.若qji≠0,則稱胞j為胞i的像胞.若令矩陣Q={qji}表示一步轉移概率矩陣,并且p(n)=[p0(n),p1(n),…,pNc(n)]T表示系統(tǒng)響應的概率分布,那么方程(6)可以寫成矩陣的形式,即：

p(n+1)=Qp(n),n=0,1,2,…

(7)

下面介紹在一個周期上運用短時高斯逼近方法計算轉移概率矩陣的方法.

1.2 短時高斯逼近解

首先考慮系統(tǒng)(1)響應過程X(t)的一階矩和二階矩:

m(t)=E[X(t)]

C(t)=E[(X-m)(X-m)T]

(8)

其中m(t)是響應過程X(t)的均值向量,C(t)是X(t)的協(xié)方差矩陣.運用It公式可得到確定一階矩m(t)和二階矩C(t)的微分方程如下:

σ(X)σT(X)]

(9)

(10)

其中h(m,C,t)和g(m,C,t)均為非線性函數(shù).

(11)

其中m(kτ)和C(kτ)是方程(10)滿足初始條件m[(k-1)τ]=x0和C[(k-1)τ]=0時的短時解,這樣能使計算結果更準確[15,16].在時間段[(k-1)τ,kτ]內,從常規(guī)胞i到胞j的轉移概率為:

(12)

(13)

p(n+1)=Q(M)Q(M-1)…Q(1)p(n),n=0,1,2,…

(14)

上述過程詳細介紹如何構造含周期激勵的系統(tǒng)在一個周期上的映射.運用方程(14)計算穩(wěn)態(tài)響應概率密度函數(shù)時的迭代終止條件為:

(15)

其中ε>0是給定的小參數(shù).

為了定量的估計方法的精度,選用直接Monte Carlo模擬方法的結果作為參考.定義如下絕對誤差積分:

(16)

其中pGCM(x,t)表示由基于短時高斯逼近的廣義胞映射方法得到的概率密度函數(shù),而pMC(x,t)是通過直接Monte Carlo模擬方法計算得到的概率密度函數(shù).

2 含指數(shù)積分型非粘性阻尼的系統(tǒng)

考慮一個周期激勵下具有指數(shù)積分型非粘性阻尼的系統(tǒng)[10],其動力學方程為:

(17)

(18)

根據(jù)Leibnitz微積分理論,方程(17)轉化為如下的三維系統(tǒng):

(19)

現(xiàn)在假設系統(tǒng)(19)還受到加性高斯白噪聲的激勵,對應的It隨機微分方程如下:

dX1=X2dt

(20)其中D表示隨機激勵的強度,B(t)是標準的Weiner過程,滿足:

E[dB(t)]=0,

(21)

如果令一階矩mi=E[Xi],二階矩vi=E[(Xi-mi)2],(i=1,2,3),cij=E[(Xi-mi)(Xj-mj)], (i,j=1,2,3,且i

(22)

矩方程(22)的初值條件為:

v1(0)=v2(0)=v3(0)=c12(0)=c13(0)=c23(0)=0

(23)

3 系統(tǒng)的穩(wěn)態(tài)響應分析

圖1 當β=0.6時,系統(tǒng)(20)的穩(wěn)態(tài)邊緣概率密度函數(shù)(a)位移X1, eIAE分別為0.0597(ζ=0.05)和0.0927(ζ=0.1),(b)速度X2, eIAE分別為0.0371(ζ=0.05)和0.0503(ζ=0.1)Fig.1 Steady-state marginal probability density functions of system (20) when β=0.6. (a)Displacement X1, where eIAE=0.0597(ζ=0.05) and eIAE=0.0927(ζ=0.1), (b)Velocity X2, where eIAE=0.0371(ζ=0.05) and eIAE=0.0503(ζ=0.1)

首先,固定松弛參數(shù)β=0.6,圖1給出了阻尼系數(shù)ζ取不同值時位移X1和速度X2的穩(wěn)態(tài)邊緣概率密度函數(shù).實線表示基于短時高斯逼近的廣義胞映射方法的結果,圓圈表示對系統(tǒng)(20)直接Monte Carlo模擬得到的結果,發(fā)現(xiàn)兩種結果基本一致.而且p(x1)的絕對誤差積分eIAE在ζ=0.05時為0.0597,在ζ=0.1時為0.0927.p(x2)的絕對誤差積分eIAE在ζ=0.05時為0.0371,在ζ=0.1時為0.0503.這些說明了短時高斯逼近和廣義胞映射方法的有效性.圖2進一步給出了由廣義胞映射方法得到的穩(wěn)態(tài)聯(lián)合概率密度函數(shù).從圖1和2中,均可以發(fā)現(xiàn)在一定范圍內概率密度函數(shù)的峰值隨著阻尼系數(shù)ζ增大而變高.

圖2 當β=0.6時,系統(tǒng)(20)的穩(wěn)態(tài)聯(lián)合概率密度函數(shù)(a)ζ=0.05, (b) ζ=0.1Fig.2 Steady-state joint probability density functions of system (20) when β=0.6(a)ζ=0.05, (b) ζ=0.1

這里需要指出,運用短時高斯逼近構造一系列的轉移概率矩陣,顯著提高了廣義胞映射方法的效率.針對該算例中轉化后的三維系統(tǒng),轉移概率矩陣的計算過程僅消耗了928s(程序運行環(huán)境為i7-6500U四核處理器筆記本電腦,編程語言為Fortran).而在相同的胞劃分情形下,利用Monte Carlo模擬方法計算轉移概率矩陣一般需要數(shù)十個小時.

然后固定阻尼系數(shù)ζ=0.15,考慮松弛參數(shù)β取不同值時系統(tǒng)的穩(wěn)態(tài)響應概率密度函數(shù).通過計算,得到系統(tǒng)位移X1和速度X2的邊緣概率密度函數(shù)如圖3所示(實線為廣義胞映射方法的結果,圓圈為直接模擬的結果),廣義胞映射方法的結果與直接模擬的結果基本吻合.而圖4則給出了對應聯(lián)合概率密度函數(shù)的廣義胞映射方法結果.觀察發(fā)現(xiàn)在一定取值范圍內松弛參數(shù)β越小穩(wěn)態(tài)概率密度函數(shù)的峰值越高.

圖3 當ζ=0.15時,系統(tǒng)(20)的穩(wěn)態(tài)邊緣概率密度函數(shù)(a)位移X1, eIAE分別為0.1236(β=0.4)和0.0986(β=0.8),(b)速度X2, eIAE分別為0.0708(β=0.4)和0.0574(β=0.8)Fig.3 Steady-state marginal probability density functions of system (20) when ζ=0.15.(a)Displacement X1, where eIAE=0.1236(β=0.4) and eIAE=0.0986(β=0.8), (b)Velocity X2, where eIAE=0.0708(β=0.4) and eIAE=0.0574(β=0.8)

圖4 當ζ=0.15時,系統(tǒng)(20)的穩(wěn)態(tài)聯(lián)合概率密度函數(shù)(a) β=0.4, (b) β=0.8Fig.4 Steady-state joint probability density functions of system (20) when ζ=0.15(a) β=0.4, (b) β=0.8

4 結論

本文運用基于短時高斯逼近的廣義胞映射方法,研究了含指數(shù)積分型非粘性阻尼和周期激勵系統(tǒng)在高斯白噪聲作用下的穩(wěn)態(tài)響應概率密度函數(shù).首先簡單介紹了方法的實施過程,然后針對一類含指數(shù)積分型非粘性阻尼的系統(tǒng),將其擴展成三維系統(tǒng),應用It公式和高斯閉包方法推導了系統(tǒng)的矩方程,計算了其穩(wěn)態(tài)響應概率密度函數(shù),并重點分析了阻尼系數(shù)和松弛參數(shù)對穩(wěn)態(tài)響應的影響.通過直接Monte Carlo模擬的結果驗證了廣義胞映射方法的有效性,并用絕對誤差積分對廣義胞映射計算結果的精度做了定量評估.從計算消耗的時間上看,短時高斯逼近方法在構造轉移概率矩陣時效率較高,因此,基于短時高斯逼近的廣義胞映射方法為研究較高維非線性系統(tǒng)的隨機響應提供了有效途徑.

1Sun J Q. Random vibration analysis of a non-linear system with dry friction damping by the short-time Gaussian cell mapping method.JournalofSoundandVibration, 1995,180(5):785～795

2Wu Y, Zhu W Q. Stochastic analysis of a pulse-type prey-predator model.PhysicalReviewE, 2008,77(4):041911

3Yue X L, Xu W, Wang L, Zhou B C. Transient and steady-state responses in a self-sustained oscillator with harmonic and bounded noise excitations.ProbabilisticEngineeringMechanics, 2012,30(5):70～76

4Yue X L, Xu W, Jia W T, Wang L. Stochastic response of a φ6 oscillator subjected to combined harmonic and Poisson white noise excitations.PhysicaA:StatisticalMechanicsanditsApplications, 2013,392(14):2988～2998

5Han Q, Xu W, Yue X L. Stochastic response analysis of noisy system with non-negative real-power restoring force by generalized cell mapping method.AppliedMathematicsandMechanics(EnglishEdition), 2015,36(3):329～336

6Li Z G, Jiang J, Hong L. Transient behaviors in noise-induced bifurcations captured by generalized cell mapping method with an evolving probabilistic vector.InternationalJournalofBifurcationandChaos, 2015,25(08):1550109

7Sun J Q, Hsu C S. The generalized cell mapping method in nonlinear random vibration based upon short-time Gaussian approximation.ASMEJournalofAppliedMechanics, 1990,57(4):1018～1025

8Han Q, Xu W, Sun J Q. Stochastic response and bifurcation of periodically driven nonlinear oscillators by the generalized cell mapping method.PhysicaA:StatisticalMechanicsanditsApplications, 2016,458:115～125

9Xu W, Rong H W, Fang T. Visco-elastic systems under both deterministic and bound random parameteric excitation.AppliedMathematicsandMechanics, 2003,24(9):1089～1099

10 Sieber J, Wagg D J, Adhikari S. On the interaction of exponential non-viscous damping with symmetric nonlinearities.JournalofSoundandVibration, 2008,314(1):1～11

11 Zhao X, Xu W, Yang Y G, Wang X Y. Stochastic responses of a viscoelastic-impact system under additive and multiplicative random excitations.CommunicationsinNonlinearScienceandNumericalSimulation, 2016,35:166～176

12 Wang D L, Xu W, Gu X D, Pei H Q. Response analysis of nonlinear vibro-impact system coupled with viscoelastic force under colored noise excitations.InternationalJournalofNon-LinearMechanics, 2016,86:55～65

13 Wu W F, Lin Y K. Cumulant-neglect closure for non-linear oscillators under random parametric and external excitations.InternationalJournalofNon-LinearMechanics, 1984,19(4):349～362

14 Sun J Q, Hsu C S. Cumulant-neglect closure method for nonlinear systems under random excitations.ASMEJournalofAppliedMechanics, 1987,54(3):649～655

15 Risken H. The Fokker Planck equation, Methods of solution and application, 2nd Ed. Berlin Heidelberg: Springer-Verlag, 1989

16 Sun J Q. Stochastic dynamics and control. Elsevier, 2006

17 Yu J S, Cai G Q, Lin Y K. A new path integration procedure based on Gauss-Legendre scheme.InternationalJournalofNon-LinearMechanics, 1997,32(4):759～768

18 Stroud A H. Numerical quadrature and solution of ordinary differential equations. New York: Springer, 1974

*The project supported by the National Natural Science Foundation of China (11532011, 11472212).

? Corresponding author E-mail:hanqun6@yahoo.com

17 March 2017,revised 18 April 2017.

STEADY-STATE STOCHASTIC RESPONSE OF A SYSTEM WITH NON-VISCOUS EXPONENTIAL DAMPING AND PERIODIC EXCITATION*

Han Qun1?Xu Wei2

(1.CollegeofScience,HuazhongAgriculturalUniversity,Wuhan430070,China) (2.DepartmentofAppliedMathematics,NorthwesternPolytechnicalUniversity,Xi′an710072,China)

The steady-state response of a system with the non-viscous exponential damping and perodic excitation is studied by the generalized cell mapping method based on short-time Gaussian approximation. The process of the method is firstly introduced, and the moment equations are derived. The steady-state response probability density functions are then presented, in which the effect of damping coefficient and relaxation parameter on the probability density functions is discussed. The validity of the generalized cell mapping method is demonstrated by the results from direct Monte Carlo simulation.

generalized cell mapping method, short-time Gaussian approximation, non-viscous exponential damping, periodic excitation, steady-state response

*國家自然科學基金資助項目(11532011, 11472212)

10.6052/1672-6553-2017-027

2017-03-17收到第1稿,2017-4-18收到修改稿.

? 通訊作者 E-mail:hanqun6@yahoo.com