多自由度粘弹性非线性随机系统的瞬态响应*

金涛 金肖玲 黄志龙†

(1.浙江省电力设计院水工结构专业, 杭州 310012)(2.浙江大学工程力学系, 杭州 310027)

多自由度粘弹性非线性随机系统的瞬态响应*

金涛1金肖玲2黄志龙2†

(1.浙江省电力设计院水工结构专业, 杭州 310012)(2.浙江大学工程力学系, 杭州 310027)

研究了高斯白噪声激励下多自由度粘弹性非线性系统的瞬态响应.首先,通过将粘弹性项对系统的作用近似地简化为对原系统阻尼部分以及刚度部分的修正,得到近似的不具粘弹性项的等效非线性随机系统. 然后,应用基于广义谐和函数的随机平均法,导出关于幅值瞬态概率密度的平均Fokker-Planck-Kolmogorov方程.该方程的解可通过多重级数式表示,基函数为幅值相关正交函数,系数为时间函数.应用Galerkin方法,关于时间的系数可由一阶线性微分方程组解得,从而得出幅值响应的瞬态概率密度、状态空间概率密度及幅值统计矩的半解析表达式.最后,以耦合的二自由度Duffing-van der Pol振子系统为例,通过与原系统数值模拟结果的比较分析验证了所提出的半解析方法的有效性,并讨论了粘弹性对系统响应的影响.

瞬态响应, 粘弹性, 多自由度非线性随机系统, 随机平均法, Galerkin方法

引言

无论是线性还是非线性动态系统,当外部载荷是高斯白噪声随机激励时,系统的响应都是扩散的马尔科夫过程,响应的转移概率密度由Fokker-Planck-Kolmogorov(FPK)方程所控制,在相应的边界条件与初始条件下求解FPK方程就可以获得系统响应的概率密度,它能够完全描述系统的响应,从而可知系统响应随时间的变化过程,完全描述其演化规律.

对于随机系统稳态响应概率密度的求解,已经发展了很多方法,如局部平衡法、平稳势法以及广义平稳势方法等等[1,2],但是能够求解瞬态响应的FPK方程则是极少数,只有一些非常特殊的一阶系统及多自由度常系数线性随机系统能够精确获得其响应的瞬态概率密度[3-5],大多数情况下只能求近似解. 至今已经有一些近似方法应用于非线性随机系统的瞬态响应研究,Atkinson基于变分法提出了FPK方程的特征值和特征函数来近似表示非线性随机系统的近似瞬态转移概率密度[6];Johnson和Scott将求解高斯白噪声下Duffing振子瞬态概率密度的问题转化为以摄动方法求特征值的问题[7];Wen通过Galerkin方法得到了单自由度非线性随机系统的近似瞬态概率密度[8].实际工程中,结构系统一般要用多自由度系统来描述,由于系统自由度的增加,高斯白噪声外激下的系统相对应的FPK方程具有更高的维数,使其求解难度更高.随机平均法,除了能保持系统的本质性质,使原来不是扩散过程的系统转化为扩散过程,还可以在一定程度上降低系统的维数,得到更低维数的FPK方程,从而降低求解难度,因此随机平均法在多自由度非线性随机响应预测中有着广泛的应用.

基于随机平均法,对单自由度非线性随机系统的近似瞬态概率密度研究已经取得一些成果,Iwan和Spanos通过等效线性化和随机平均法给出了高斯白噪声激励下非线性系统关于幅值瞬态概率密度的FPK方程[9],进而以摄动法求得以特征值和特征向量近似表示的幅值瞬态概率密度;Spanos等人将随机平均法应用于高斯白噪声激励下的单自由度非线性随机系统,将该系统的幅值瞬态解近似地表示为Rayleigh分布与适当选取正交基函数的线性组合,组合系数是随时间变化的并可通过Galerkin方法得到[10].而对于多自由度非线性随机系统,主要通过等效线性化方法或Monte Carlo方法对其进行瞬态响应统计量的研究[11-15],对系统响应的瞬态概率密度的研究还比较少,文献[16]基于广义谐和函数的随机平均法将Spanos等[10]提出的关于幅值瞬态响应半解析方法成功地推广于多自由度的非线性随机系统,文献[17]进一步应用于多自由度时滞非线性随机系统的瞬态响应分析.

考虑到实际工程中的材料有存在粘弹性,本文研究多自由度粘弹性非线性随机系统的瞬态响应.首先将含粘弹性项的系统近似等效为不具粘弹性项的系统.然后应用基于广义谐和函数的随机平均法,导出幅值概率密度的FPK方程,应用文献[16,17]的方法,得到瞬态响应的半解析解.最后以耦合的二自由度Duffing-van der Pol振子系统为典型算例验证求解过程.

1 无粘弹性的等效系统

由于复合材料、聚合物以及合金的广泛应用,这些材料的粘弹性特性对于结构系统特性带来较大的影响,同时在随机外激或参激激励下,这些结构系统在一定程度上可由多自由度粘弹性非线性随机系统来描述,设其动力学控制方程为：

(i,j=1,…,n;k=1,…,m)

(1)

当ε=0时,系统(1)可简化为多自由度的哈密顿系统,其控制方程为：

(2)

该系统每个自由度的首次积分为：

(3)

上式中Ui(xi)=∫xi0g(ui)dui为势能函数.假设系统(2)的解为原点平衡点附近的周期解,则它们可表示为[18]：

xi(t)=aicosθi(t)

θi(t)=φi(t)+φi(t)

(i=1,...,n)

(4)

上式中ai为瞬时幅值,θi为瞬时相位,vi为瞬时频率且可表示为：

vi(ai,θi)=dφi/dt=

(5)

则平均频率可以表示为：

(6)

通过哈密顿函数和势能之间的关系Hi=Ui(ai),平均频率也可写为：

(7)

当ε是小量的情况下,Xi(t-τ)中的τ不是特别大的时候,可以用如下表达式近似[19]：

(8)

利用h(t)函数的指数衰减特性,结合式(8),系统(1)中的松弛算子可以近似简化为：

(9)

上式中包含时间衰减项e-(t/λil)随时间以指数快速衰减,因此该项在近似分析时可以被忽略,粘弹性项对系统(1)的影响简化为两部分,其一是对系统刚度的修正,另一部分是对系统阻尼的修正,从而原系统(1)可近似等效为：

(i,j=1,…,n;k=1,…,m)

(10)

该近似等效系统中的阻尼项与刚度项分别为：

(11)

式(11)中,δij为Kronecker delta符号下同.

2 系统的近似瞬态响应

2.1 随机平均法

系统(1)的响应可通过对系统(10)的研究来近似,由于系统(10)中每个振子具有弱阻尼、弱激励的性质,同时各个振子之间刚度是不耦合的,因此其广义速度和广义位移可以通过广义谐和函数的形式来表示：

Xi(t)=Ai(t)cosΘi(t)

(i=1,…,n)

(12)

类似于式(4)与式(5),瞬时相位Θi和瞬时频率Λi可表示为：

Θi(t)=Φi(t)+Γi(t)

Λi(Ai,Θi)=dΦi/dt=

(13)其中Ai(t),Θi(t),Φi(t)以及Γi(t)为随机过程,将式(12)代入近似系统(10)中并利用式(12)的协调条件,可得关于幅值Ai和相位Γi的随机微分方程：

dAi/dt=εF1i(A,Γ)dt+ε1/2G1ik(A,Γ)ξk(t)

dΓi/dt=εF2i(A,Γ)dt+ε1/2G2ik(A,Γ)ξk(t)

(i=1,…,n;k=1,…,m)

(14)

其中A=[A1,…,An]T为幅值向量,Γ=[Γ1,…,Γn]T为相位向量,式中各项参数为：

c″ij(A,Θ)AjΛj(Aj,Θj)sinΘj

c″ij(A,Θ)AjΛj(Aj,Θj)sinΘj

(15)

dAi=mi(A)dt+σik(A)dBk(t)

(i=1,…,n;k=1,…,m)

(16)

上式中的平均漂移与平均扩散系数为：

bij(A)=σik(A)σjk(A)=ε〈2DkkG1ikG1jk〉Θ

(17)

(18)

为简便起见,假设原系统(1)开始处于静止状态,则式(18)的初始条件为：

(19)

2.2 基于Laguerre多项式的近似解

一般很难得到方程(18)的精确解析解,这里用近似方法求解.将FPK方程(18)的解,即幅值概率密度p(A,t),近似地表示为[16]：

p(A,t)=

(20)

上式可理解为相应线性系统的解与非线性使解偏离线性系统之解此两部分之和,其中sr1,…,rn(t)为待定的随时间变化组合系数,ηiri和Riri(Ai)分别为系统(1)中第i个非耦合线性振子系统对应的FPK方程的特征值和特征函数[16]：

(i=1,…,n;ri=0,1,…)

(21)

上式中Lri(·)为第ri阶Laguerre多项式,2Dlili代表子系统的外激励强度.特征函数Riri(Ai)具有如下一些基本性质：

Ri(r+1)(Ai)=

(22)

将近似表达式(20)代入FPK方程(18),可得残差：

(23)

利用Galerkin方法,可令误差R在一定条件下取得最小值来得到线性组合系数sr1,…,rn(t),这里取如下的条件：

(k1,…,kn=0,1,…)

(24)

结合上式与式(23),并利用Laguerre多项式和特征函数的基本性质(22),可得关于线性组合系数sr1,…,rn(t)的一阶线性常微分方程组：

(k1,…kn=0,1,…)

(25)

将近似表达式(20)代入式(19),并利用Laguerre多项式的性质可得微分方程组(25)的初始条件：

sr1,…,rn(0)=0 (r1,…,rn=0,1,2,…)

(26)

通过数值计算求得sr1,…,rn(t)并代入式(20)可得近似幅值概率密度p(A,t),实际计算时对近似解(20)中的级数作适当截断,即取ri=0,1,…Ni(i=1,…,n),使其满足一定的精度要求.则系统幅值响应的稳态概率密度求解简化为求以下线性代数方程组相应矩阵的特征值和特征向量问题：

ki=0,…,Ni(i=0,1,…,n)

(27)

幅值稳态概率密度可由上式中的未知量Sr1,…,rn来表示：

(28)

(29)

其相应的广义位移和广义速度瞬态联合概率密度：

(31)

Il,r=0,(r>l/2)

(32)

其中l为偶数.

3 算例

为了验证求解方法的可行性及精度,考虑高斯白噪声激励下两个Duffing-Van der Pol耦合振子系统,其动力学控制方程为：

(33)

(34)

修正后的阻尼及刚度系数为：

(35)

根据式(5)及式(6),平均频率ωa1可通过下式以数值方法求得：

ωai(Ai)=

(36)

其中Elliptick(·)为第一类完全椭圆积分.应用2.1节中的随机平均法,可得关于幅值Ai的平均It随机微分方程,形如式(16),相应的平均FPK方程形如式(18),根据式(17)导出其平均漂移与平均扩散系数为：

b12=b21=0

(37)

利用2.2节提出的近似解表示方法,近似幅值瞬态概率密度可以表示为：

R1i(A1)R2j(A2)

(38)

其中相应线性系统的特征值和特征函数为：

(k=1,2)

(39)

应用Galerkin方法可得如下一阶微分方程组：

(k=0,…,N1,l=0,…,N2)

(40)

利用特征函数的正交性,上式中各项满足如下关系：

fm(Am)=

gm(Am)=

Iij=(2j+1)δij-jδi,j-1-(j+1)δi,j+1

Bi(·)=Li(·)/i!

(m=1,2;i=0,…N1;j=0,…N2)

(41)

在初始条件(26)下求解方程组(40),将结果代入(38)可得近似幅值瞬态概率p(A1,A2,t),从而可以得到每个振子的近似幅值瞬态概率密度p(A1,t)与p(A2,t),其具体形式为：

(42)

系统响应幅值的近似稳态概率密度由式(27)与(28)导出.由式(30)可得每个振子的广义位移和广义速度联合瞬态概率密度：

(k=1,2)

(43)

基于幅值瞬态概率密度(42)以及幅值统计矩(31),幅值的一阶和二阶统计矩具有如下形式：

(44)

图1 系统第一个振子幅值瞬态概率密度p(A1,t); ps(A1)为幅值稳态概率密度;—本文方法的结果;●,○,▼表示对原系统(33)模拟结果Fig. 1 The transient probability densities of amplitude response for the first oscillator, where is the stationary probability density, ‘-’ indicates the results from the proposed procedure, ●,○,▼ presents the simulation results of the original system (33)

图2 系统第一个振子幅值的一阶与二阶统计矩E[A1(t)]和随时间的演化过程;—本文方法的结果;●对原系统(33)模拟结果Fig. 2 The time evolution of the first- and second- statistical moments‘—’ indicates the results from the proposed procedure; first oscillator,●,○,▼ present the simulation results of the original system (33)

图3 系统第一个振子广义速度和广义位移联合概率密度1,t)的演化过程为稳态概率密度Fig. 3 The evolution of the transient joint probabilitydensities is the stationary probability density

图4 不同粘弹性参数E11下系统第一个振子幅值一阶与二阶统计矩E[A1(t)]和随时间的演化过程;—本文方法的结果;●,○,▼,▽对原系统(33)模拟结果Fig. 4 The time evolution of the first- and second-statistical moments oscillator with different values of E11, where ‘—’ indicates the results from the proposed procedure, ●,○,▼,▽present the simulation results of the original system (33)

图5 不同粘弹性参数λ11下系统第一个振子幅值一阶与二阶统计矩 E[A1(t)]和随时间的演化过程;—本文方法的结果; ●,○,▼对原系统(33)模拟结果Fig. 5 The time evolutions of the first- and second-statistical moments oscillator with different values of λ11, where ‘—’ indicates the results from the proposed procedure, ●,○,▼present the simulation results of the original system (33)

为了反映粘弹性对系统响应的影响,图4和图5给出了不同粘弹性参数下一阶和二阶统计矩随时间的演化.从图4可以看出,当粘弹性参数E11逐渐增大时,第一个振子的幅值响应逐渐减小,这可以从修正阻尼和修正刚度的表达式(35)中来看出,随着E11的增大,第一个振子的修正阻尼系数变大而修正刚度系数变小,从而其响应减小.还可以看出随着时间的增大,粘弹性参数对系统响应的影响越来越明显.从图5可以看出,当粘弹性参数λ11逐渐增大时,第一个振子的幅值响应也是逐渐减小的.

需要指出的是,由式(9)对粘弹性项进行简化过程中,e-(t/λil)中λil的取值量级为10-1,当时间t较大时,该项很快衰减,因此它只对初期的很短时间具有较大影响,本文考虑瞬态时间t大于10的结果,其受到的影响不大,这可从所得到的计算结果看出.

4 结论

本文研究了多自由度粘弹性非线性随机系统瞬态响应.提出了多自由度粘弹性非线性随机系统的随机平均法,将粘弹性对系统的作用近似地简化为对系统阻尼以及刚度的修正,从而将原系统近似等效为无粘弹性项的多自由度非线性随机系统,应用基于广义谐和函数的随机平均导出关于幅值响应瞬态概率的平均FPK方程.幅值瞬态解以半解析的方式近似地表示多重级数式表示,基函数为幅值相关正交函数,所选取的基函数以Laguerre多项式为主要组成部分,并且是与原系统相对应线性系统的FPK方程的特征函数,组合系数是随时间变化的,并可通过Galerkin方法得到.以非线性阻尼耦合的Duffing-Van der Pol振子系统为例,分析幅值瞬态概率密度,广义位移与广义速度联合瞬态概率密度及幅值的一二阶统计矩并讨论粘弹性对响应的影响,本文得到的近似半解析解与Monte Carlo数值模拟结果吻合很好,说明本文的近似方法可有效地应用于多自由度粘弹性非线性随机系统响应的研究.

1 朱位秋. 随机振动. 北京：科学出版社,1992 (Zhu W Q. Random vibration. Beijing: Science Press, 1992(in Chinese))

2 朱位秋. 非线性随机动力学与控制: Hamilton 理论体系框架. 北京: 科学出版社,2003 (Zhu W Q. Nonlinear stochastic dynamics and control : Hamiltonian formulation. Beijing: Science Press, 3003(in Chinese))

3 Caughey T K. Nonlinear theory of random vibrations.AdvancesinAppliedMechanics, 1971,11:209～253

4 Caughey T K, Dienes J K. Analysis of a nonlinear first‐order system with a white noise input.JournalofAppliedPhysics, 1961,32(11):2476～2479

5 Gardiner C W. Handbook of stochastic methods for physics, chemistry and the natural sciences. Berlin: Springer, 1983

6 Atkinson J D. Eigenfunction expansions for randomly excited non-Linear systems.JournalofSoundandVibration, 1973,30(2):153～172

7 Johnson J P, Scott R A. Extension of eigenfunction-expansion solutions of a Fokker-Planck equation-II. second order System.InternationalJournalofNon-LinearMechanics, 1980,15(1):41～56

8 Wen Y K. Approximate method for nonlinear random vibration.JournaloftheEngineeringMechanicsDivision, 1975,101(4):389～401

9 Iwan W D, Spanos P D. Response envelope statistics for nonlinear oscillators with random excitation.ASMEJournalofAppliedMechanics, 1978,45(1):170～174

10 Spanos P D, Sofi A, Di Paola M. Nonstationary response envelope probability densities of nonlinear oscillators.ASMEJournalofAppliedMechanics, 2007,74(2):315～324

11 ASME Micaletti R C, et al. A solution method for linear and geometrically nonlinear MDOF systems with random properties subject to random excitation.ProbabilisticEngineeringMechanics, 1998,13(2):85～95

12 Zhang L, Zu J W, Zheng Z. The stochastic Newmark algorithm for random analysis of multi-degree-of-freedom nonlinear systems.Computers&Structures, 1999,70(5):557～568

13 Ohtori Y, Spencer Jr B F. Semi-implicit integration algorithm for stochastic analysis of multi-degree-of-freedom structures.JournalofEngineeringMechanics, 2002, 128(6): 635～643

14 Smyth A W, Masri S F. Nonstationary response of nonlinear systems using equivalent linearization with a compact analytical form of the excitation process.ProbabilisticEngineeringMechanics, 2002,17(1):97～108

15 Saha N, Roy D. The girsanov linearization method for stochastically driven nonlinear oscillators.ASMEJournalofAppliedMechanics, 2007,74(5):885～897

16 金肖玲. 多自由度强非线性随机系统的响应与稳定性研究[博士学位论文]. 杭州：浙江大学,2009 (Jin X L. Response and stability of multi-degree-of-freedom strongly nonlinear stochastic systems[PhD Thesis]. Hangzhou: Zhejiang University, 2009 (in Chinese))

17 金涛, 黄志龙. 调制白噪声激励下时滞非线性系统的瞬态响应. 中国科学: 物理学力学 天文学, 2013,43:1～8 (Jin T, Huang Z L. Transient probability densities of nonlinear system with time delay subject to modulated white noise excitation.ScienceChinaPhysics,Mechanics&Astronomy, 2013,43:1～8 (in Chinese))

18 Xu Z, Cheung Y K. Averaging method using generalized harmonic gunctions for strongly non-Linear oscillators.JournalofSoundandVibration, 1994,174(4):563～576

19 Liu Z H, Zhu W Q. Stochastic averaging of quasi-integrable Hamiltonian systems with delayed feedback control.JournalofSoundandVibration, 2007,299(1):178～195

20 Stratonovich R L. Topics in the theory of random noise: Vol.1: general theory of random processes: nonlinear transformations of signals and noise. New York: Gordon and Breach, 1963

21 Khasminskii R Z. A limit theorem for the solutions of differential equations with random right-hand sides.TheoryofProbability&ItsApplications, 1966, 11(3):390～406

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

† Corresponding author E-mail: zlhuang@zju.edu.cn

3 April 2017,revised 18 April 2017.

TRANSIENT RESPONSE OF NONLINEAR MULTI-DEGREE-OF-FREEDOM STOCHASTIC SYSTEM WITH VISCOELASTICITY*

Jin Tao1Jin Xiaoling2Huang Zhilong2†

(1.ZhejiangElectricPowerDesignInstitute,Hangzhou310012,China) (2.DepartmentofEngineeringMechanics,ZhejiangUniversity,Hangzhou310027,China)

The transient response of a nonlinear multi-degree-of-freedom system with viscoelasticity subjected to Gaussian white noise excitation is investigated. Firstly, the effect of the viscoelasticity on the system is approximated by the modified damping and stiffness. The original system is replaced by a system without viscoelasticity. Then, the stochastic averaging method based on generalized harmonic functions is adopted to derive the averaged Fokker-Planck-Kolmogorov equation of amplitude transient joint probability density for each oscillator. This equation is solved by expressing the probability density as multiple series in terms of a set of properly state-dependent orthogonal basis functions with time-dependent coefficients. According to Galerkin method, the time-dependent coefficients can be solved from a set of first-order linear differential equations. Finally, the semi-analytical formulae of the transient probability density as well as the transient probability of the state response and the statistical moments for the amplitude response is obtained. To illustrate the proposed procedure, coupled two-degree-of- freedom Duffing-van der Pol oscillators with viscoelasticity subjected to Gaussian white noise excitation is investigated as an example. The effect of viscoelasticity on the system response is initially discussed. Moreover, comparison with the simulation results of the original system indicates that the proposed procedure is accuracy and efficacy.

transient response, viscoelasticity, nonlinear multi-degree-of-freedom stochastic system, stochastic averaging, Galerkin method

*国家自然科学基金资助项目(11532011, 11672262, 11621062)

10.6052/1672-6553-2017-025

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

† 通讯作者 E-mail: zlhuang@zju.edu.cn