基于Hertz接触的单自由度碰振系统的随机响应近似闭合解

祝海生,陈林聪,孙建桥,赵珧冰

(1.华侨大学 土木工程学院，福建 厦门 361021；2.加州大学Merced分校 工程学院，美国加利福尼亚州 95343；3.天津大学 机械工程学院，天津 300072)

碰撞振动系统普遍存在于工程领域中。由于碰撞振动系统的强非线性，使得关于碰撞振动系统的一些动力学性能研究变得极其复杂，如碰振系统中的随机响应、分叉以及混沌等。实际的工程领域中的一些碰振现象亦体现了碰撞振动系统产生的影响，如桥梁梁端的碰撞振动导致的梁端开裂；海洋中轮船与冰块的碰撞振动造成的船体损伤；列车车轮与铁轨之间的碰撞所造成的一系列的材料、能量消耗。因此，基于碰撞振动现象的普遍性以及在各工程领域产生的影响，对于碰撞振动的研究具有十分重要的实际意义。

碰撞振动是一类非常复杂的动力学过程。根据碰撞过程中所表现的一系列动力学和物理学特征，例如碰撞时物体的相对速度、碰撞接触时间、碰撞接触面的形状变化以及碰撞区域的弹塑性变形等。碰撞振动现象的描述一般主要有两种模型，分别是瞬时冲击模型和Hertz接触理论及其修正模型。目前，基于瞬时冲击模型，Dimentberg[1-2]就随机激励下线性碰撞振动系统的随机响应分析分别做了大量的研究。王亮等[3]定义了一种随机响应度量研究了一类二自由度碰撞振动系统在随机白噪声激励下的随机响应，通过与数值模拟比较讨论了随机噪声对于系统的影响。田海勇等[4]通过建立一类单自由度含间隙的碰撞振动系统的动力学模型，讨论了随机干扰对碰撞振动系统的影响。伍新等[5]考虑碰撞振动系统的Poincare映射的隐式特点研究了一类三自由度含间隙双面碰撞振动系统Poincare映射的叉式分岔的反控制问题。Gu等[6]提出了一种预测碰撞系统响应的随机平均法，利用所提出的随机平均法考察了高斯白噪声激励下的单自由度碰撞振动系统，得到了关于系统能量的平稳概率密度函数。类似的，Liu等[7]讨论了碰撞系统在有色噪声激励下的平稳概率响应，详细地讨论了有色噪声和恢复系数对碰撞系统响应的影响。Li等[8]则获得了相关高斯白噪声激励下Duffing-van der Pol碰撞振动系统的平稳概率密度函数，考察了对不同参数引起的随机分岔问题。徐伟等[9]先借助非光滑变换和狄拉克函数，随后得到了Duffing-Rayleigh 碰撞振动系统的对应等效非线性系统，最后采用奇异摄动法分析了等效非线性系统的随机P-分岔。Zhu[10-11]分别将高斯白噪声和泊松白噪声激励下Duffing碰撞振动系统的稳态响应概率密度近似为指数多项式形式。与蒙特卡罗模拟结果进行对比，该指数多项式形式的近似解具有较高的精度，但不能反映系统局部的一些非光滑本质特征。Chen等[12]应用迭代加权残值法获得了高斯白噪声激励下单自由度碰撞振动系统的稳态响应闭合解。

Hertz接触理论模型既能够很好的描述在碰撞过程中的物体的局部变形，同时也能够较好的反映碰撞过程中的接触力的变化。在20世纪90年代，Jing等[13-14]基于Hertz接触理论模型率先获得了高斯白噪声外激情形时的单自由度系统的精确平稳解。后来，Huang等[15]利用Hertz理论模型将约束模型转化为于非线性弹簧，考察了二自由度碰撞振动系统的随机稳态响应，并在一定条件下得到了这类系统的精确平稳解。Xu等[16-17]联合利用等效非线性方法和随机平均法将Hertz模型进行改进分析了单自由度非弹性碰撞振动系统的响应。另外，徐明等[18]采用修正Hertz模型基于随机平均法分析了在高斯白噪声激励下非弹性碰撞振动系统的首次穿越问题，得到了系统的条件可靠性函数和相应的条件概率密度函数。

特别需要指出的是，当前的碰撞随机振动研究还具有一些局限性。如随机平均法须在系统恢复系数接近于1以及在弱阻尼弱激励等情况下才能应用；等效线性化法能够较为准确的得到系统的均方速度和均方位移，然而其应用范围却往往局限于高斯统计的情形，并且在参激情况下常被认为是不充足或不合适的；摄动法只适用于小参数、弱非线性系统条件下，对于非线性较强的情况下的误差非常大，甚至其结果根本不正确。总之，关于随机碰撞振动的研究，特别是在随机响应闭合解方面仍需投入大量的研究。

本文应用随机振动研究的最新成果——迭代加权残值法，求解了高斯白噪声激励下基于Hertz接触的单自由度强非线性碰振系统的随机响应近似闭合解。首先，利用概率环流和概率势流的概念，构造系统对应简化FPK方程解的近似表达式；然后，运用加权残值法获得系统平稳概率密度函数的近似闭合解；最后，应用迭代的办法，提高近似闭合解的精度。研究结果表明，本文的迭代加权残值法采用分段函数形式表达，具有较高的精度，可以有效的解决非光滑问题，同时还清楚的表征了碰振系统的非光滑本质特征等等。作为算例，分别研究Duffing碰振系统和干摩擦碰振系统，来验证本文方法的有效性。

1 问题概述

考虑一个受高斯白噪声激励的单自由度Hertz碰振系统，其图示和运动方程如图1所示。

图1 碰撞振动模型Fig.1 Vibro-impact system

(1)

式中，g1(X,Y)表示线性或非线性阻尼系数，g2(X)表示系统非碰撞回复力，f(X)表示碰撞回复力，其表达式为

(2)

式中，系统平衡点与两侧弹性壁的相互作用规律满足Hertz接触定律，δr和δl是碰撞振动系统平衡点距两端碰撞壁的距离，Br和Bl是常数，其与弹性壁的材料及几何形状有关；hi(X,Y)表示线性或非线性激励幅值，Wi(t)表示高斯白噪声，其相关函数为E[Wi(t)Wj(t+τ)] =2Dijδ(τ),i,j=1,2,…,l。

与系统(1)相对应的简化FPK方程如下

(3)

式中，p=p(x,y)为式(1)的平稳概率密度函数；漂移系数(m1,m2)和扩散系数(b22)如下

m1=y

m2=g1(x,y)y+g2(x)+f(x)-

(4)

根据详细平衡法[19-20]，式(3)可分为以下两个部分

(5)

和

(6)

其中式(5)描述了概率环流的平衡，表示保守系统的能量守恒。式(6)则表明了概率势流的平衡，表示为系统能量输入与耗散之间的平衡。可见，详细平衡法具有明确的物理意义。若系统存在精确平稳解，式(5)和式(6)必同时被满足。然而，在实际工程中，该情形常难以实现。因此用详细平衡法求解FPK方程的精确平稳解的能力和范围是有限的。

迭代加权残值法[21]是近期作者提出来的一种求解FPK方程平稳解的有效方法。当系统存在精确解时，近似闭合解收敛于精确解。当系统不存在精确解，获得具有较高精度的近似闭合解。本文将推广该方法求解简化FPK方程(3)。

2 迭代加权残值法

迭代加权残值法主要由加权残值法和迭代步骤组成。在具体的应用过程中，先结合概率势流和概率环流概念，构造简化FPK方程 (3)的近似表达式；再用加权残值法获得近似闭合解；最后引入迭代技术来提高近似闭合解的精度。

2.1 加权残值法

首先，假定方程(3)的平稳解是如下形式

(7)

式中,C0是归一化常数，φ(x,y)关于状态变量的多项式，按照文献[21]可由以下形式三部分组成

k2ψ2(x,y)

(8)

式中,cij,k1,k2是待求系数，ψ1和ψ2是概率势流和概率环流，分别满足式(5)与(6)，表达式如下

lnb22(x,y)

(9)

φ(x,y)→∞ asA→∞,α>-1

ifφ(x,y)∝AαasA→0

(10)

式中,x=Asinθ，y=Acosθ，其中第二个条件保证了在原点的可积性。

(11)

(12)

式中,Ml(x,y)(l=1,2,…,N+2,0

Mk(x,y)=pm(x,y)xiyj,k=1,2,…,N

MN+1(x,y)=pm(x,y)ψ1(x,y)

MN+2(x,y)=pm(x,y)ψ2(x,y)

(13)

式中,pm可为由等效线性化获得的高斯概率分布函数或随机平均法获得的概率分布函数，也可构造为如下形式的函数

pm(x,y)=exp[-ψ1(x,y)-ψ2(x,y)]

(14)

2.2 迭代过程

(15)

为了评价近似解的精度，将近似解与参考解pR(x,y)之间的均方误差的平方根定义为

(16)

式中，pR=pR(x,y)可以是蒙特卡罗模拟结果，也可以是精确解析解。

需指出的是，迭代方向的选择也是至关重要的。若迭代方向选取的不好，可能会获得一个不收敛的结果。这一现象在一些强非线性系统当中表现的尤为突出。因此，Chen等[23]提出了一种渐近迭代的概念，具体步骤是先应用加权残值法计算弱非线性系统或弱激励情形时的近似平稳概率密度函数，然后以此构造权函数，再在非线性参数空间里，逐步搜索强非线性系统的近似平稳概率密度函数。

(17)

接下来，我们将分别考察两个算例，说明迭代加权残值法在Hertz碰撞随机振动系统研究中的适用性。

3 算 例

3.1 算例1

考虑在高斯白噪声激励下的Duffing碰撞随机振动系统，运动方程如下

(18)

式中，β1>0，Wi(t)表示激励强度为2Di的独立的高斯白噪声。

对应FPK方程的漂移系数和扩散系数如下

(19)

由此可得出概率势流和概率环流为以下形式

(20)

下面分别考察了对称情形与非对称情形。

对称情形：设δr=δl=0.5，β1=0.1，D1=0.1，α1=1，α3=0.5，Br=Bl=10，经2次迭代之后所得结果如下

p(x,y)=0.228 791 97×exp[-0.942 889 34x2+

0.007 802 02xy+0.180 384 93y2-

0.075 531 07x4+0.0008 465 5x3y-

0.196 890 08x2y2-0.002 027 02xy3-

0.378 371 96×

(21)

根据式(16)，式(21)的精度为0.013 3。

图2给出了系统(18)在Br=Bl=10,δr=δl=0.5对称情形时的平稳响应概率密度函数。其中图2(a)与(b)分别表示其联合概率密度函数的理论近似闭合解(21)与相应的蒙特卡罗模拟结果(样本数为108)；图2(c)与(d)分别表示关于位移的边缘概率密度函数与关于速度的边缘概率密度函数，其中实线表示解析结果，符号(○)表示蒙特卡罗模拟结果。如图所示，迭代加权残值法计算的结果与蒙特卡罗模拟结果非常相近，精度令人满意；系统平稳响应概率密度函数是单峰态，且顶部光滑；初步表明了本文方法在Duffing系统对称情形中是有效的。

(a)

(b)

(c)

(d)图2 系统(18)在Br=Bl=10，δr=δl=0.5情形时的平稳概率密度函数Fig.2 The stationary PDF of system (18)in the case of Br=Bl=10 and δr=δl=0.5

非对称情形：为了进一步考察迭代加权残值的有效性和精度，现以平稳响应概率密度函数(21)构造权函数，应用迭代加权残值法分别考虑不同弹性壁刚度Br(Bl)与不同碰撞壁间距参数δr(δl)等非对称情形下的平稳响应概率密度函数的分布。

首先，针对不同弹性壁刚度Br(Bl)情形，现选取Br=100,Bl=10。经3次迭代得到精度为0.020 5的结果如下

p(x,y)=0.262 593 24×exp[0.101 677 44x-

0.002 757 12y-0.914 625 26x2+

0.008 178 59xy+0.055 220 25y2-

0.736 233 78x3+0.000 004 81x2y-

0.007 625 82xy2+0.006 756 1y3-

0.494 431 02x4-0.003 171 1x3y-

0.167 376 57x2y2-0.002 123 55xy3-

0.497 403 68×

(22)

类似地，图3给出了系统(18)在Br=100，Bl=10情形时的平稳响应概率密度函数。由图3可知，迭代加权残值法计算结果依然能够吻合的很好；系统平稳响应概率密度函数整体呈光滑“竹笋”状，关于位移的边缘概率密度具有显著的非对称性；迭代加权残值法进一步得到了有效验证。另外，由图2和图3可以得出，系统关于速度的边缘概率密度函数函数基本相同。

现结合对称情形，并增加计算Br=20，Bl=10这一情况(受篇幅所限，不再赘述具体表达式，下同)，考虑系统在一侧弹性壁刚度不同的情形下对于位移的边缘概率密度函数的影响情况。由下图4可知,随着弹性壁刚度Br(Bl)的增大，关于位移的边缘概率密度函数峰值会变大，系统稳态响应会有所减小。

(a)

(b)

(c)

(d)图3 系统(18)在Br=100,Bl=10情形时的平稳响应概率密度函数，其他参数同图2Fig.3 The stationary PDF of system (18)in the case of Br=100，Bl=10.The other parameters are the same as those in Fig.2

图4 系统(18)分别在Br=Bl=10、Br=20,Bl=10以及Br=100,Bl=10时关于位移边缘概率密度函数，其中实线表示理论解析解，符号(o,△,□ )表示蒙特卡罗模拟结果,其他参数同图2Fig.3 The marginal PDF of displacement of system (18)in the cases of Br=Bl=10,Br=20,Bl=10 or Br=100,Bl=10,respectively.Solid line represents the analytical solution,symbols(o,△,□ )represent the Monte Carlo simulation data.The other parameters are the same as those in Fig.2

然后固定参数Br(Bl)，考察不同碰撞壁间距参数δr(δl)情形，现选取δr=1，δl=0.5,Br=Bl=10，经4次迭代得到精度为0.025 1的结果如下

p(x,y)=0.195 716 30×exp[-0.058 343 74x+

0.002 590 17y-0.677 146 06x2+

0.010 489 24xy+0.200 581 69y2+

0.181 499 36x3+0.000 701 53x2y+

0.012 265 64xy2-0.000 749 63y3-

0.123 542 70x4+0.000 511 15x3y-

0.330 673 92×

0.167 784 26x2y2-0.002 800 01xy3-

(23)

由图5可知，在不同碰撞壁间距参数下应用迭代加权残值法计算的结果与样本数为108的蒙特卡罗模拟结果吻合的效果也能够令人满意，关于位移的边缘概率密度数也具有较明显的非对称性。类似于图4，图6讨论了系统平衡点到碰撞壁一侧间距增大时关于位移边缘概率密度函数的分布情况。由图6可知，随着系统平衡点到碰撞壁一侧间距增大时，系统的“单峰”单侧坡度会有所变缓，系统的稳态响应会增大。

3.2 算例2

算例1研究表明，迭代加权残值法在光滑Hertz碰振系统研究中可取得较好的效果。下面将考察一类非光滑系统Hertz碰撞系统，来进一步检验迭代加权残值法的适用性。现以干摩擦碰振系统为例，其运动方程如下

(a)

(b)

(c)

(d)图5 系统(18)在δr=1,δl=0.5时的平稳概率密度函数，其他参数同图2Fig.5 The stationary PDF of system (18)in the case of δr=1,δl=0.5.The other parameters are the same as those in Fig.2

图6 系统(18)在δr=δl=0.5,δr=1,δl=0.5与δr=2.5,δl=0.5时关于位移边缘概率密度函数。其中实线表示理论解析解，(o,△,□ )表示蒙特卡罗模拟结果，其他参数同图2Fig.6 The marginal PDF of displacement of system (18)in the cases of δr=δl=0.5,δr=1,δl=0.5 or δr=2.5,δl=0.5，respectively.The solid line represents the analytical solution,symbols (o,△,□ )represent the Monte Carlo simulation data.The other parameters are the same as those in Fig.2

(24)

β0,β1,β3≥0，sgn(Y)为符号函数，W(t)表示激励强度为2D1的高斯白噪声。系统的平稳FPK方程形如式(3)，对应的漂移系数和扩散系数如下

b22=2D1

(25)

相应概率势流和概率环流如下形式

(26)

同样，针对干摩擦系统亦考虑对称情形与非对称情形。

对称情形：取参数：δr=δl=0.5,Br=Bl=5,β0=β1=β3=0.1,D1=0.1。应用迭代加权残值法经过3次迭代计算得到精度为0.018 1的结果如下

p(x,y)=0.280 673 82×exp[0.049 948 83x2-

0.009 358 99xy+0.463 128 39y2-

0.191 118 95x4+0.007 645 05x3y+

0.018 861 58x2y2+0.012 918 00xy3-

1.087 645 1×

0.262 478 15y4-1.082 915 1|y|]

(27)

图7为系统(24)在Br=Bl=5对称情形时的平稳响应概率密度函数结果。图7(a)与(b)分别表示基于式(27)的联合稳态概率密度函数与相应的蒙特卡罗模拟结果(样本数为108)。图7(c)与(d)分别表示关于位移的边缘概率密度函数与关于速度的边缘概率密度函数。由图7可知，理论解析结果与蒙特卡罗模拟结果非常相近；系统平稳响应概率密度函数是单峰，顶部呈“刀锋”状，充分体现出了干摩擦系统的非光滑特征。

(a)

(b)

(c)

(d)图7 系统(24)在Br=Bl=5时的平稳响应概率密度函数Fig.7 The stationary PDF of system (24)in the symmetry case of Br=Bl=5

接着，以式(27)构造权函数，经3次迭代，求得了Br=Bl=10情形下精度为0.020 3的平稳响应联合概率密度函数，其表达式如下

p(x,y)=0.319 755 06×exp[0.145 562 56x2-

0.014 088 440xy+0.464 592 29y2-

0.509 706 45x4+0.016 608 98x3y+

0.032 655 97x2y2+0.014 236 69xy3-

0.992 554 01×

0.272 600 96y4-1.097 960 4|y|]

(28)

根据式(28)以及蒙特卡罗模拟结果得出了系统在Br=Bl=10情形时的平稳响应概率密度函数，如图8所示。由图8可知，与Br=Bl=5情形相比，Br=Bl=10情形时系统“刀锋”的开口宽度有所减小；系统的干摩擦特征仍获得了比较好的体现。此外，比较图7和图8(以及下图10和图11)也可以得出，系统在所研究的情形下关于速度的边缘概率密度函数也几乎相同。

为了考察不同弹性壁刚度Br(Bl)对系统响应的影响，我们也考察了Br=Bl=20这一情况。图9给出了不同弹性壁刚度情形下关于位移的边缘概率密度函数对比图。由图9知，当系统弹性壁刚度逐渐增大时，系统的“刀锋”开口宽度是逐渐减小的，系统的稳态响应也在减小；另外，初步得出系统的“刀锋”宽度与系统的弹性壁刚度大小是无关的。

由以上可知，应用迭代加权残值法在干摩擦系统对称情形中，当弹性壁刚度逐步增大时，仍然能够保持较高的精度，初步表明了本文方法在干摩擦碰撞振动系统对称情形中的有效性和适用性。

(a)

(b)

(c)

(d)图8 系统(24)在Br=Bl=10时的平稳概率密度函数，其他参数与图7相同Fig.8 The stationary PDF of system (24)in the case of Br=Bl=10.The other parameters are the same as those in Fig.7

图9 系统(24)分别在Br=Bl=5、Br=Bl=10与Br=Bl=20时关于位移的边缘概率密度函数。其中实线表示理论解析解，(o,△,□ )表示蒙特卡罗模拟结果Fig.9 The marginal PDF of displacement of system (24)in the cases of Br=Bl=5,Br=Bl=10 and Br=Bl=20,respectively.The solid line represents the analytical solution,symbols(o,△,□ )represent the Monte Carlo simulation data

非对称情形：现固定参数Br=Bl=10，考察不同碰撞壁间距参数δr(δl)情形。先取δr=1.5,δl=0.5,β0=β1=β3=0.1,D1=0.1。以式(28)构造权函数，经3次迭代之后得出精度为0.015 3的结果。表达式如下：

p(x,y)=0.192 588 83×exp[-0.000 033 18x+

0.000 439 14y+0.001 236 61x2-

0.001 225 29xy+0.493 140 42y2+

0.000 024 08x3+0.000 068 03x2y-

0.000 062 91xy2-0.000 213 25y3-

0.000 912 73x2y2+0.002 382 99xy3-

1.256 498 2×

0.000 074 00x4+0.000 030 95x3y-

0.250 424 29y4-1.013 793 0|y|]

(29)

图10进一步给出了δr=1.5,δl=0.5非对称情形时的平稳响应概率密度函数。由图10可知，关于位移的边缘概率密度函数亦表现出明显的非对称性，系统平稳响应概率密度函数的峰值区间基本维持在1.5～-0.5；另外，迭代加权残值法计算的结果与样本数为108的蒙特卡罗模拟结果吻合的也非常好，干摩擦的非光滑性也得到了充分体现。

现保持δl不变，取δr=2。以式(29)构造权函数，经2次迭代获得精度为0.014 0的平稳响应联合概率密度函数闭合解如下

(a)

(b)

(c)

(d)图10 系统(24)在δr=1.5,δl=0.5时的平稳概率密度函数Fig.10 The stationary PDF of system (24)in the case of δr=1.5,δl=0.5

p(x,y)=0.163 760 83×exp[-0.010 293 78x-

0.006 261 64y+0.027 450 11x2-

0.001 538 97xy+0.492 225 00y2+

0.002 910 81x3+0.000 761 82x2y+

0.010 068 95xy2+0.005 455 80y3-

0.005 365 30x4-0.001 042 16x3y-

0.022 606 36x2y2+0.003 372 81xy3-

1.354 765 5×

0.250 615 05y4-1.018 964 6|y|]

(30)

图11给出了系统在δr=2.0,δl=0.5情形时的平稳响应概率密度函数。此时由于δr的增大，系统平稳响应概率密度函数的峰值区间也基本保持在2.0～-0.5；结合图10和图11可以发现，系统的“刀锋”开口宽度随着碰撞壁间距(δr)的增大有所增大。图12给出了δr=1.0,δl=0.5、δr=1.5,δl=0.5与δr=2，δl=0.5情形时关于位移的边缘概率密度函数对比图，考察了不同碰撞壁间距系数δr(δl)对系统响应的影响。由图12可知，所举三种情况下的平稳响应边缘概率密度函数峰值区间基本都维持在x=δr和x=δl之间，再结合图9和图12可以得出，系统的“刀锋”的区间宽度只与δr(δl)有关，与Br(Bl)无关；另外，边缘概率密度函数p1(x)峰值分别在0.42、0.35和0.3左右，比较三种情况得知，当系统平衡点到碰撞壁一侧间距增大时，关于位移的边缘概率密度的“刀锋”峰值会有所降低，同时“刀锋”开口宽度会有所增大；研究表明应用迭代加权残值法在干摩擦碰撞振动系统的研究中也是有效的。

(a)

(b)

(c)

(d)图11 系统(24)在δr=2,δl=0.5时的平稳概率密度函数，其他参数同图10Fig.11 The stationary PDF of system (24)in the case of δr=2,δl=0.5.The other parameters are the same as those in Fig.10

图12 系统(24)分别在δr=1.0,δl=0.5、δr=1.5,δl=0.5与δr=2,δl=0.5时关于位移的边缘概率密度函数。其中实线表示理论解析解，符号(o,△,□ )表示蒙特卡罗模拟结果。其他参数同图10Fig.12 The marginal PDF of displacement of system (24)in the cases with δr=1.0,δl=0.5、δr=1.5,δl=0.5 or δr=2,δl=0.5,respectively.The solid line represents the analytical solution,symbols(o,△,()represent the Monte Carlo simulation data.The other parameters are the same as those in Fig.10

4 结 论

本文采用迭代加权残值法研究了随机激励下基于Hertz接触的单自由度强非线性碰撞振动系统平稳响应问题。研究方法主要有三个步骤：首先构造支配系统平稳响应概率密度函数的简化FPK方程的近似表达式；然后，应用加权残值法获得系统平稳概率密度函数的近似闭合解；最后应用迭代技术，来提高近似闭合解的精度。作为算例，分别研究了Duffing碰振系统以及干摩擦碰振系统。研究表明本文的方法在光滑系统(Duffing碰振系统)以及非光滑系统(干摩擦碰振系统)研究中均是有效和适用的，同时还完全表征了系统的非线性动力学特征，如干摩擦系统顶部的“刀锋”非光滑特征。