含指数积分型阻尼和周期激励系统的稳态随机响应分析*

韩群 徐伟

(1.华中农业大学理学院,武汉 430070) (2.西北工业大学应用数学系,西安 710072)

含指数积分型阻尼和周期激励系统的稳态随机响应分析*

韩群1†徐伟2

(1.华中农业大学理学院,武汉 430070) (2.西北工业大学应用数学系,西安 710072)

运用基于短时高斯逼近的广义胞映射方法,研究了含指数积分型非粘性阻尼和周期激励系统在高斯白噪声作用下的稳态响应.首先介绍了方法的实施过程,并推导了系统的矩方程.然后给出了系统的稳态概率密度函数,分析了阻尼系数和松弛参数对稳态响应的影响,并通过直接Monte Carlo模拟的结果验证了广义胞映射方法的有效性.

广义胞映射, 短时高斯逼近, 指数积分型阻尼, 周期激励, 稳态响应

引言

广义胞映射方法已被广泛用来研究非线性动力系统的随机响应分析[1-6].传统广义胞映射在计算状态胞之间一步转移概率时,通常借助Monte Carlo数值模拟,对于含周期激励的系统,在考虑一步转移概率时,转移时间长度直接取为系统的周期,虽然这样在形式上显得很简单,但是运算过程却相当费时.因此,该方法在含周期激励高维系统中的应用显得非常困难.Sun和Hsu[7]提出一种短时高斯逼近方法构造FPK方程的短时解,从而计算了高斯白噪声激励下典型自治系统的一步转移概率矩阵,该方法显著提高广义胞映射方法的效率,在应用中体现出了较大的优势.最近,Han等[8]将基于短时高斯逼近的广义胞映射方法推广,用于研究周期和高斯白噪声激励下非线性系统的响应概率密度函数,并研究了一类光滑非连续(SD)振子的瞬态和稳态响应行为.

指数积分型非粘性阻尼是一种复杂的阻尼模型[9,10],与经典粘性阻尼不同,它不仅与当前的瞬时速度有关,还依赖于速度的历史状态.含指数积分型非粘性阻尼系统的随机响应研究中,已有一些成果[11,12],但是这些都没有考虑周期激励的作用.在处理含这类阻尼的系统时,通常将阻尼看作是一个变量,由此将系统扩展成三维系统,若考虑周期和高斯白噪声共同激励下该系统的随机响应分析,必然会对传统的广义胞映射方法的应用提出更大的挑战.本文基于文献[8]中广义胞映射方法,讨论了含指数积分型非粘性阻尼和周期激励系统在高斯白噪声作用下的稳态响应分析.首先简单介绍含周期激励系统随机响应分析的广义胞映射方法,以及短时高斯逼近在一个周期上的实施策略.然后针对一类含指数积分型非粘性阻尼的系统,先将其扩成三维系统,给出了系统的矩方程.最后计算了系统的稳态响应概率密度函数,并讨论了阻尼有关参数对稳态响应的影响.

1 基于短时高斯逼近的广义胞映射求解过程

考虑含周期和高斯白噪声激励的N维非线性随机动力系统,其对应的It随机微分方程为：

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

(1)

其中X(t)=[X1(t),X2(t),…,XN(t)]T∈RN(上标T表示转置)是一个N维Markov随机过程,代表系统的响应.f(X,t)是一个N×1的漂移向量函数,其中包含周期激励,即满足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)的形式导数即为L维的标准高斯白噪声.

1.1 广义胞映射方法

n=0,1,2,…

(3)

在广义胞映射方法的实现过程中,连续的状态空间RN被转化成离散的胞状态空间ZN.考虑一个感兴趣的有界区域Ω,并将之划分为Nc个常规胞,它们的尺寸大小相同,假设为h1×h2×…×hN,其中hi是每个胞在xi方向上的长度.所有的常规胞依次编号为1到Nc, 区域Ω以外的所有部分被看作是一个陷胞,编号为0. 如果pi(n)表示t=nT时系统的响应位于第i号胞的概率,而qji表示以pi(0)=1为初始条件,系统的响应在T时刻位于第j号胞的转移概率,那么有：

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

(4)

对于一个常规胞i,qji可由如下公式计算得到：

(5)

(6)

该方程确定了一个有限Markov链的演化.若qji≠0,则称胞j为胞i的像胞.若令矩阵Q={qji}表示一步转移概率矩阵,并且p(n)=[p0(n),p1(n),…,pNc(n)]T表示系统响应的概率分布,那么方程(6)可以写成矩阵的形式,即：

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

(7)

下面介绍在一个周期上运用短时高斯逼近方法计算转移概率矩阵的方法.

1.2 短时高斯逼近解

首先考虑系统(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)的协方差矩阵.运用It公式可得到确定一阶矩m(t)和二阶矩C(t)的微分方程如下:

σ(X)σT(X)]

(9)

(10)

其中h(m,C,t)和g(m,C,t)均为非线性函数.

(11)

其中m(kτ)和C(kτ)是方程(10)满足初始条件m[(k-1)τ]=x0和C[(k-1)τ]=0时的短时解,这样能使计算结果更准确[15,16].在时间段[(k-1)τ,kτ]内,从常规胞i到胞j的转移概率为:

(12)

(13)

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

(14)

上述过程详细介绍如何构造含周期激励的系统在一个周期上的映射.运用方程(14)计算稳态响应概率密度函数时的迭代终止条件为:

(15)

其中ε>0是给定的小参数.

为了定量的估计方法的精度,选用直接Monte Carlo模拟方法的结果作为参考.定义如下绝对误差积分:

(16)

其中pGCM(x,t)表示由基于短时高斯逼近的广义胞映射方法得到的概率密度函数,而pMC(x,t)是通过直接Monte Carlo模拟方法计算得到的概率密度函数.

2 含指数积分型非粘性阻尼的系统

考虑一个周期激励下具有指数积分型非粘性阻尼的系统[10],其动力学方程为:

(17)

(18)

根据Leibnitz微积分理论,方程(17)转化为如下的三维系统:

(19)

现在假设系统(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 系统的稳态响应分析

图1 当β=0.6时,系统(20)的稳态边缘概率密度函数(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)

首先,固定松弛参数β=0.6,图1给出了阻尼系数ζ取不同值时位移X1和速度X2的稳态边缘概率密度函数.实线表示基于短时高斯逼近的广义胞映射方法的结果,圆圈表示对系统(20)直接Monte Carlo模拟得到的结果,发现两种结果基本一致.而且p(x1)的绝对误差积分eIAE在ζ=0.05时为0.0597,在ζ=0.1时为0.0927.p(x2)的绝对误差积分eIAE在ζ=0.05时为0.0371,在ζ=0.1时为0.0503.这些说明了短时高斯逼近和广义胞映射方法的有效性.图2进一步给出了由广义胞映射方法得到的稳态联合概率密度函数.从图1和2中,均可以发现在一定范围内概率密度函数的峰值随着阻尼系数ζ增大而变高.

图2 当β=0.6时,系统(20)的稳态联合概率密度函数(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

这里需要指出,运用短时高斯逼近构造一系列的转移概率矩阵,显著提高了广义胞映射方法的效率.针对该算例中转化后的三维系统,转移概率矩阵的计算过程仅消耗了928s(程序运行环境为i7-6500U四核处理器笔记本电脑,编程语言为Fortran).而在相同的胞划分情形下,利用Monte Carlo模拟方法计算转移概率矩阵一般需要数十个小时.

然后固定阻尼系数ζ=0.15,考虑松弛参数β取不同值时系统的稳态响应概率密度函数.通过计算,得到系统位移X1和速度X2的边缘概率密度函数如图3所示(实线为广义胞映射方法的结果,圆圈为直接模拟的结果),广义胞映射方法的结果与直接模拟的结果基本吻合.而图4则给出了对应联合概率密度函数的广义胞映射方法结果.观察发现在一定取值范围内松弛参数β越小稳态概率密度函数的峰值越高.

图3 当ζ=0.15时,系统(20)的稳态边缘概率密度函数(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时,系统(20)的稳态联合概率密度函数(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 结论

本文运用基于短时高斯逼近的广义胞映射方法,研究了含指数积分型非粘性阻尼和周期激励系统在高斯白噪声作用下的稳态响应概率密度函数.首先简单介绍了方法的实施过程,然后针对一类含指数积分型非粘性阻尼的系统,将其扩展成三维系统,应用It公式和高斯闭包方法推导了系统的矩方程,计算了其稳态响应概率密度函数,并重点分析了阻尼系数和松弛参数对稳态响应的影响.通过直接Monte Carlo模拟的结果验证了广义胞映射方法的有效性,并用绝对误差积分对广义胞映射计算结果的精度做了定量评估.从计算消耗的时间上看,短时高斯逼近方法在构造转移概率矩阵时效率较高,因此,基于短时高斯逼近的广义胞映射方法为研究较高维非线性系统的随机响应提供了有效途径.

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