随机激励下滞迟系统的稳态响应闭合解1)

刘 俊 陈林聪,2) 孙建桥

∗(华侨大学土木工程学院，福建厦门361021)†(加州大学Merced分校工程学院，美国加利福尼亚州95343)

-动力学与控制

随机激励下滞迟系统的稳态响应闭合解1)

刘 俊∗陈林聪∗,2)孙建桥†

∗(华侨大学土木工程学院，福建厦门361021)†(加州大学Merced分校工程学院，美国加利福尼亚州95343)

滞迟系统属于一类典型的强非线性系统，滞迟力不仅取决于系统的瞬时变形，还与变形历程有关.虽然滞迟系统的随机振动问题已被广泛研究，但至今尚未得到滞迟系统随机响应概率密度函数的精确闭合解.本文运用迭代加权残值法获得了高斯白噪声激励下Bouc-Wen滞迟系统稳态响应概率密度函数的近似闭合解.首先，运用等效线性化法求出系统的稳态高斯概率密度函数；然后以此构造权函数，应用加权残值法求得了系统指数多项式形式的非高斯概率密度函数；最后引入迭代的过程，逐步优化权函数，提高计算所得结果的精度.以随机地震激励下钢纤维陶粒混凝土结构的稳态响应作为算例，其中Bouc-Wen模型的参数是基于拟静力学试验数据，并应用最小二乘法辨识获得.与Monte Carlo模拟结果相比，等效线性化法得到的结果精度较差；由加权残值法得到的结果能够表现出非线性特征，但其精度依然无法令人满意；采用迭代加权残值法得到的近似闭合解与Monte Carlo模拟的结果吻合非常好；对于较强随机激励情形，采用渐进迭代加权残值法具有较高的求解效率，所获得的理论解析解具有较高的精度.结果表明,所获得的近似闭合解不仅对于土木工程领域具有重要的实际应用价值，而且还可作为检验其他非线性系统随机响应预测方法的精度的标准.

滞迟系统，稳态响应，闭合解，迭代加权残值法，钢纤维陶粒混凝土结构

引言

滞迟现象如弹塑性[13]、铁电性[4]、形状记忆合金材料[5]等常出现在科学研究和工程实际的不同领域，强烈震动载荷下的结构系统通常表现出滞迟现象[67].学术界提出了许多模型描述这些滞迟关系，如双线性模型[89]，Ramberg-Osgood模型[10]，Bouc-Wen模型[67]，Ozdemir模型[11]，Masing模型[12]，Duhem 模型[13]，Preisach模型[14]等，其中Bouc-Wen模型是较为通用的一种.

滞迟力不仅与系统当前的状态有关，而且还与系统过去的状态有关.因此，滞迟动力学系统属于一类典型的强非线性系统.近年来，滞迟系统的随机振动问题已被广泛研究，但至今尚未获得精确的闭合解，因此，学术界提出了各种近似解法.如，Iwan[15]采用基于克雷洛夫-包哥留波夫假设的等效线性化法研究了双线性滞迟系统的随机响应.Roberts[16-17]分别采用标准随机平均法及能量包线随机平均法研究了双线性滞迟系统的随机响应.Bouc[7]釆用FPK方程法研究了滞迟系统的随机响应.Wen[7,18]采用线性化方法研究了滞迟系统的随机响应.Zhu等[19]采用能量包线随机平均法研究分布弹塑性元件为滞迟恢复力模型的滞迟系统的随机响应.Lin等[20]采用能量包线随机平均法研究了Bouc-Wen滞迟系统的随机响应.Ying等[21]采用能量包线随机平均法研究了Duhem滞迟系统的随机响应.对具有非局部记忆特性的Preisach滞迟系统，Mayergoyz和Korman[22]研究了随机输入下 Preisach系统的平均输出，Ni等[23]基于方差和非局部记忆滞迟本构模型的切换概率分析，近似得到了Preisach滞迟系统的非线性随机动力响应的二阶统计量.Spanos等[24]和Wang等[25]分别采用标准随机平均和能量包线随机平均法研究了 Preisach滞迟系统的随机响应.最近，Jin等[26]运用随机平均法和伽辽金法获得了随机地震激励下滞迟系统的近似瞬态响应.

然而，上述方法都存在一些不足，如随机平均法仅限于能量耗散很小且弱激励情形；等效线性化法能得到较准确的均方速度和均方位移，但是这种方法只局限于高斯统计的情形，在参激情况下常被认为是不充足或不合适的；非高斯闭合法在尾部通常会获得负值概率，特别是参激情形.因此，还需进一步开展滞迟非线性系统的随机响应预测研究.

Er[2729]提出了一种指数多项式闭合法求解稳态FPK方程，但由于采用了高斯概率密度函数构造权函数，因此指数多项式闭合解法的应用范围受到了很大限制.Di Paola和Sof[30]改进了指数多项式闭合法，提出采用一种简单有效的迭代过程以提高近似解的精确性，但该方法没有从根本上解决指数闭合解的局限性问题.最近，文献[31]提出了一种求解FPK方程稳态解的新方法——迭代加权残值法，该方法的核心思想是逐步优化权函数，目前已被成功应用于求解多种复杂单自由度强非线性系统[3132].

本文将迭代加权残值法进一步应用于Bouc-Wen滞迟系统，构造其稳态概率密度函数.首先，将稳态FPK方程的解设为指数多项式，然后借助加权残值法确定假设解中的待定系数；引入迭代方法，逐步优化权函数，提高加权残值法的精度，得到系统稳态响应的近似闭合解.作为应用，本文获得了高斯白噪声激励下钢纤维陶粒混凝土框架结构的近似稳态响应概率密度函数闭合解.研究表明，所获得的闭合解与Monte Carlo模拟结果吻合得较好.

1 迭代加权残值法

迭代加权残值法主要由两个步骤组成.首先，应用加权残值法获得系统的近似稳态概率密度响应；然后，引入迭代过程，逐步优化权函数，提高加权残值法的精度，最终得到具有较高精度的近似闭合稳态解.

1.1 加权残值法

考虑外激高斯白噪声激励下的 Bouc-Wen系统，其运动方程如下

式中，X,Y与Z分别表示系统位移、速度与滞迟力；ξ为黏滞阻尼率；α∈(0,1)是屈服前后的刚度比；λ，β，γ与n为滞回环参数，β与γ控制滞回环的形状，λ控制滞回环的幅值，n控制滞回环曲线的光滑性；W(t)是强度为2Dδ(t)的高斯白噪声.

与系统(1)相应的稳态FPK方程为

式中，p=p(x,y,z).目前，式(2)尚未获得精确解析解.

构造如下形式的近似解

其中，C0为归一化常数，φ(x,y,z,cijk)为关于状态变量的n阶多项式，可表示为

其中，cijk为待求系数.¯p(x,y,z)的存在条件为

其中，x=Rsinθcosφ,y=Rsinθsinφ,z=Rcosθ.将方程(3)代入式(2)中，得残差

由于¯p(x,y,z)只是p(x,y,z)的近似值.因此，残差r(x,y,z,cijk)通常不为零.根据加权残值法，引入一组权函数 Mijk(x,y,z)，使其与残差r(x,y,z,cijk)的乘积在域内的积分为零，即

权函数可取为[27-29]

其中,pm(x,y,z)是由等效线性化法或随机平均法等方法获得的系统概率密度函数.值得注意的是，在数值计算时方程(7)中的积分域通常为有限积分域，但通常难以确定.为此，运用Monte Carlo模拟法粗略估计积分域Ωs.方程(7)可近似表示为

其中，0＜i+j+k≤l,l为解的阶数.将式(10)进行数值积分，即可得到一组关于cijk的二次非线性代数方程组.采用牛顿法求解该方程组，可获得系统的近似稳态响应概率密度.

1.2 迭代过程

滞迟系统属于强非线性系统，仅使用一次加权残值法可能得不到具有足够精度的解.因此，采用文献[30-31]中的迭代技术来提高解的精度.令(k)为迭代k次后得到的近似稳态概率密度函数，取代式(9)中的pm，再应用加权残值法求出下一个近似稳态概率密度函数¯p(k+1).重复此步骤直到满足以下收敛条件

其中，ε0为预设误差，N1,N2与N3为常数，分别表示状态空间离散的数目，Δx,Δy和Δz表示离散积分步长.

其中，pR(x,y,z)表示Monte Carlo模拟得到的稳态概率密度函数.

需要指出的是,当pm与真实解相差较大时，上述迭代过程可能不会收敛，特别是当系统为强非线性的情况.文献[31]提出了一种渐近迭代的方案，即先应用加权残值法获得系统弱非线性情形时或弱阻尼弱激励情形时的近似稳态概率密度函数，然后以此构造权函数，再在非线性参数空间渐近迭代.例如，针对滞迟系统(1)，为了获得激励强度D=0.2情形时的近似稳态响应概率密度函数,首先采用迭代加权残值法获得弱随机激励D=0.1情形时的近似稳态概率函数pm1；然后将pm1代入式(9)，以此构造新的权函数，再次利用迭代加权残值法，求得随机激励D=0.15情形时的近似稳态概率函数 pm2；最后将pm2代入式(9)，构造新的权函数，利用迭代加权残值法求解D=0.2时的概率密度函数.这种渐近迭代的方法可有效避免选取的初始值 pm与真实解相差较大时导致的不收敛情况，提高了求解的效率.

2 算例

考虑如图1所示的水平随机地震作用下单自由度钢纤维陶粒混凝土框架滞迟系统，其平衡条件为

其中，x(t)是集中质量的相对位移；xt(t)是集中质量的总位移；g(x,)是恢复力.

图1 地震地面运动下单自由度钢钎维陶粒混凝土框架结构Fig.1 Steel fibe reinforced ceramsite concrete(SFRCC)single degree of freedom frame subjected to earthquake ground motion

由图1中可得

式中，xg和分别表示地面的位移和加速度,将方程(15)代入方程(13)可得

其中，F(t)为图2所示导致水平地面加速度的等效载荷.该框架结构被简化为如图2所示的等效系统.

图2 等效系统Fig.2 Equivalent system

将方程(17)无量纲化，并把等效载荷理想化为独立的高斯白噪声.该系统的动力学方程化为

其中，x为无量纲化的位移,ξ为黏滞阻尼率,α是屈服前后的刚度比,k1为常数,W1(t)是强度为2D1的高斯白噪声,z为恢复力中依赖于时间的滞迟力，由以下公式表示

本文基于钢纤维陶粒混凝土柱受水平低周反复加载试验所得的数据，经最小二乘法辨识获得滞迟系统的参数：n= 1,λ=1.1037,β=0.202和γ=0.3147.理论曲线和试验曲线对比如图3所示，二者吻合较好.系统的其他参数为：ξ=0.025,α=0.2,D=0.1,ε0=0.001.

图3 钢纤维陶粒混凝土柱完全塑性阶段试验曲线和理论曲线(实线表示理论结果，虚线表示试验结果)Fig.3 Theoretical curve and practical curve of SFRCC column in fully plastic stage(the solid line denotes the theoretical curve and the dotted line denotes the practical curve)

最终可得如下形式的运动方程

本文首先运用等效线性化法获得了系统的近似稳态响应概率密度函数，表达式如下

然后以式(21)构造权函数，运用加权残值法获得了新的近似稳态响应概率密度函数，即

根据误差公式(11)，式(21)与式(22)的误差分别为0.163和0.249.显然，式(22)的精度低于式(21).因此，为了进一步提高精度，现以解(22)构造权函数，经过1次迭代，获得了误差为0.035的解

图4给出了关于p(x,z)的稳态边缘概率密度函数.由图4和图5可知，等效线性化的结果与Monte Carlo模拟结果相差甚远；加权残值法的计算结果精度尚未令人满意，但已经表现出系统的一些非线性特征；迭代加权残值法得到结果与Monte Carlo模拟的结果吻合较好.

另外，以式(23)构造权函数，采用渐进迭代加权残值法，经过3次迭代获得了较强随机激励D=0.2情形时系统的近似稳态概率密度函数闭合解

图4 D=0.1情形时关于p(x,z)的稳态边缘概率密度函数Fig.4 The steady-state marginal probability density function p(x,z)in case of D=0.1

图5 D=0.1情形时的稳态边缘概率密度函数p1(x)和p2(z)(◦,∗表示Monte Carlo模拟数据)Fig.5 The steady-state marginal probability density function p1(x)and p2(z)in case of D=0.1(◦,∗ represent the Monte Carlo simulation data)

式(24)的精度为0.040.图6给出D=0.2情形时的稳态边缘概率密度函数.图6(a)与6(b)表示关于p(x,z)的稳态边缘概率密度函数.图6(c)与图6(d)分别表示关于p1(x)及p2(z)的稳态边缘概率密度及均方差.符号(◦,∗)表示样本数为40000的Monte Carlo模拟结果.由图6可知，理论解析解与模拟结果在较强随机激励情形时仍吻合得较好.

图6 D=0.2情形时稳态边缘概率密度函数(◦,∗表示Monte Carlo模拟结果)Fig.6 The steady-state marginal probability density function in case of D=0.2(◦,∗ represent results from Monte Carlo simulation)

3 结论

本文首次获得Bouc-Wen滞迟系统的稳态响应概率密度函数的近似闭合解.首先，利用等效线性化法获得的高斯概率密度函数构造权函数，然后应用加权残值法获得系统指数多项式形式的非高斯概率密度函数，最后引入迭代的办法提高了加权残值法计算所得结果的精度.对较强随机激励情形，以弱随机激励情形时的稳态响应概率密度函数构造权函，再在随机激励的参数空间内渐近迭代.作为算例，本文研究了随机激励下钢纤维陶粒混凝土结构的稳态响应，其中Bouc-Wen系统的参数是基于拟静力学试验由最小二乘法辨识获得.对比Monte Carlo模拟结果，等效线性化的结果精度较差，加权残值法的结果虽表现出了系统的非线性特征但精度尚未令人满意，迭代加权残值法得到结果与Monte Carlo模拟的结果吻合得较好.对于较强随机激励情形，基于渐近迭代加权残值法具有较高的效率，且所得结果也具有较好的精度.本文所获得的闭合解，不仅对于土木工程领域具有重要的实际应用价值，还可作为一个标准，用来检验其他非线性系统随机响应预测方法的精度.

1万征，姚仰平，孟达.复杂加载下混凝土的弹塑性本构模型.力学学报,2016,48(5)：1159-1171(Wan Zheng,Yao Yangping,MengDa.An elastoplastic constitutive model of concrete of complicated load.Chinese Journal of Theoretical and Applied Mechanics,2016,48(5)：1159-1171(in Chinese))

2李龙彪.纤维增强陶瓷基复合材料疲劳迟滞回线模型研究.力学学报,2014,46(5)：710-729(Li Longbiao.Investigation on fatigue hysteresis loops models of fibre-reinforce ceramic-matrix composites.Chinese Journal of Theoretical and Applied Mechanics,2014,46(5)：710-729(in Chinese))

3郭洪宝，贾普荣，王波等.基于迟滞行为的2D-SiC/SiC复合材料组份力学性能分析.力学学报，2015，47(2)：260-269(Guo Hongbao,Jia Purong,Wang Bo et al.Study on constituent properties of a 2D-SiC/SiC composite by hysteresis measurements.Chinese Journal of Theoretical and Applied Mechanics,2015,47(2)：260-269(in Chinese))

4万强，陈常青，沈亚鹏.铁电陶瓷PZT53复杂力电耦合行为的实验研究.力学学报，2005,37(4)：413-420(Wan Qiang,Chen Changqing,Shen Yapeng.An experimental investigation into the complex electromechanical behavior of PZT53.Chinese Journal of Theoretical and Applied Mechanics,2005,37(4)：413-420(in Chinese))

5 Mayergoyz I.Mathematical models of hysteresis.IEEE Transactions on Magnetic,1986,22(5)：603-608

6 Wen YK.Methods of random vibrationfor inelastic structures.Applied Mechanics Review,1989,42(2)：39-52

7 Bouc R.Forced vibration of mechanical systems with hysteresis//Proceedings of the Fourth Conference on Non-Linear Oscillation,Prague,Czechoslovakia,1967

8 Caughey T.Random excitation of a system with bilinear hysteresis.ASME Journal of Applied Mechanics,1960,27(4)：649-652

9 Suzuki Y,Minai R.Application of stochastic di ff erential equations to seismic reliability analysis of hysteretic structures.Probabilistic Engineering Mechanics,1988,3(1)：43-52

10 Jennings PC.Periodic response of a general yielding structure.Journal of Engineering Mechanics-ASCE Division,1964,90(2)：131-166

11 Bhatti M,Pister K.A dual criteria approach for optimal-design for earthquake-resistant structural systems.Earthquake Engineering Structure Design,1981,9(6)：557-572

12 Beck JL,Jayakumar P.Class of Masing models for plastic hysteresis in structures//Proceedings of 14th ASCE Structures Congress,1996：1083-1090

13 Visintin A.Di ff erential Models of HysteresisⅢ.Springer Science&Business Media,2013

14 Mayergoyz I.Mathematical models of hysteresis.IEEE Transactions on Magnetic,1986,22(5)：603-608

15 Iwan WD.A distributed-element model for hysteresis and its steadystate dynamic response.ASME Journal of Applied Mechanics,1966,33：893

16 Roberts JB.The response of an oscillator with bilinear hysteresis to stationary random excitation.ASME Journal of Applied Mechanics,1978,45(4)：923-928

17 Roberts JB.Application of averaging methods to randomly excited hysteretic systems//Nonlinear Stochastic Dynamic Engineering Systems,Berlin,Springer,1988

18 Wen YK.Equivalent linearization for hysteretic systems under random excitation.ASME Journal of Applied Mechanics,1980,47(1)：150-154

19 Zhu WQ,Lei Y.Stochastic averaging of energy envelope of bilinear hysteretic systems//Nonlinear Stochastic Dynamic Engineering Systems，Berlin：Springer,1988

20 LinYK,Cai GQ.Probabilistic Structural Dynamics：AdvancedTheory and Applications.New York：McGraw-Hill,1995

21 ZG Ying，WQ Zhu，Ni YQ,et al.Stochastic averaging of Duhem hysteretic systems.Journal of Sound and Vibration,2002,254(1)：91-104

22 Mayergoyz ID,Korman CE.The Preisach model with stochastic input as a model for aftere ff ect.Journal of Applied Physics,1994,75(10)：5478-5480

23 NiYQ,YingZG,KoJM.RandomresponseanalysisofPreisachhysteretic systems with symmetric weight distribution.ASME Journal of Applied Mechanics,2002,69(2)：171-178

24 Spanos PD,Cacciola P,Muscolino G.Stochastic averaging of Preisach hysteretic systems.ASCE Journal of Engineering Mechanics,2004,130(11)：1257-1267

25 Wang Y,Ying ZQ,Zhu WQ.Stochastic averaging of energy envelope of Preisach hysteretic systems.Journal of Sound and Vibration,2009,321(3-5)：976-993

26 Jin XL,Wang Y,Huang ZL.Analysis and control for transient responses of seismic-excited hysteretic structures.Soil Dynamics and Earthquake Engineering,2015,73(6)：58-65

27 Er GK.Exponential closure method for some randomly excited non-linear systems.International Journal of Nonlinear Mechanics,2000,35(1)：69-78

28 Er GK.A consistent method for the solution to reduced FPK equation in statistical mechanics.Physica A,1999,262(1)：118-128

29 Er GK.An improved closure method for analysis of nonlinear stochastic systems.Nonlinear Dynamics,1998,17(3)：285-297

30 Di Paola M,Sof A.Approximate solution of the Fokker-Planck-Kolmogorov equation.Probabilistic Engineering Mechanics,2002,17(4)：369-384

31 Chen LC,Sun JQ.The closed-form solution of the reduced Fokker-Planck-Kolmogorov equation for nonlinear systems.Communications in Nonlinear Science&Numerical Simulation,2016,41(12)：1-10

32 Chen LC,Sun JQ,The closed-form steady-state probability density function of van der pol oscillator under random excitations.Journal of Applied Nonlinear Dynamics,2016,5(4)：495-502

THE CLOSED-FORM SOLUTION OF STEADY STATE RESPONSE OF HYSTERETIC SYSTEM UNDER STOCHASTIC EXCITATION1)

Liu Jun∗Chen Lincong∗,2)Sun Jian-Qiao†∗(College of civil Engineering，Huaqiao University，Xiamen 361021，Fujian，China)†(School of Engineering University of California Merced，CA 95343，USA)

The hysteretic system is one of the typical strongly nonlinear systems.Hysteretic force depends not only on the instantaneous deformation but also on the past history of deformation.In the last few decades,random vibration of hysteretic system has been studied extensively,but no closed-form solution of random response of hysteretic systems is available so far.In this paper,the newly developed nonlinear random vibration scheme called iterative method of weighted residuals is explored to obtain the closed-form solution of steady-state probability density function(PDF)of the Bouc-Wen hysteretic system under Gaussian white noise excitation.First,a Gaussian PDF is obtained with equivalent linearization technique,which is used as a weighting function.Then,the method of weighted residuals is utilized to determine the non-Gaussian PDF of exponential polynomial type.Finally,an iterative procedure is introduced to improvethe accuracy of the solutions obtained from the method of weighted residuals.As an illustrative example,the steadystate stochastic response of the steel fibe reinforced ceramsite concrete column under random excitation is studied,in which the hysteretic parameters associated with Bouc-Wen hysteretic model is identifie from the pseudo-static test by using the method of least square.Compared to the Monte Carlo results,the accuracy of results obtained from equivalent linearization method is poor.The results obtained from weight residue method can show the nonlinearity of Bouc-Wen systems,but its accuracy is still unsatisfactory.The iterative method of weight residuals can lead to results with higher accuracy.In the case of stronger random excitation,the progressive iterative method of weighted residuals has high efficiency.The obtained solutions agree well with the Monte Carlo simulation data.The proposed closed-form solution of PDF of Bouc-Wen hysteretic system not only is significan to the civil engineering,but also can be a benchmark to examine the accuracy of solutions obtained by other methods.

hysteretic system,steady-state response,closed-form solution,iterative method of weight residuals，steel fibe reinforced ceramsite concrete structure

O324

：A

10.6052/0459-1879-17-003

2017–01–11 收稿，2017–04–11 录用,2017–04–11 网络版发表.

1)国家自然科学基金(11172197，11332008，11572215，11672111，51608211)、福建省自然科学基金(2013J05080)、福建省高校杰出青年科研人才培育计划和华侨大学优秀青年科技创新人才(ZQN-YX307)资助项目.

2)陈林聪，副教授，主要研究方向：非线性随机振动与控制.E-mail：lincongchen@hqu.edu.cn

刘俊,陈林聪,孙建桥.随机激励下滞迟系统的稳态响应闭合解.力学学报,2017,49(3)：685-692

Liu Jun,Chen Lincong,Sun Jian-Qiao.The closed-form solution of steady state response of hysteretic system under stochastic excita-tion.Chinese Journal of Theoretical and Applied Mechanics,2017,49(3)：685-692