随机振动范文

随机振动范文（精选8篇）

随机振动 第1篇

卫星光通信与微波通信相比具有通信容量大、终端体积小和保密性好等特点,是未来卫星通信技术的一个很有潜力的发展方向。随着卫星光通信技术研究的不断深入,卫星平台随机振动的影响得到广泛的重视。有研究表明,随机振动会增大光通信的误码率[1,2],影响通信质量。

为抵偿卫星平台微振动造成的影响,振动补偿技术的研究显得尤为重要。有研究者经过深入调研后指出,导致卫星平台振动的原因繁多,机理复杂,但是对光通信影响较大的主要是频率低于100 Hz、振幅100µrad以内的振动[3]。为抵消振动给通信终端带来的种种不良影响,有研究者提出了采用面阵CCD作为探测器件测量振动,采用快速偏转镜作为执行器件的振动补偿方案,实验取得了显著的效果[3,4]。

但是,在以往的研究工作中,将卫星平台的振动简化为正弦振动进行理论分析和实验研究,还不能真实地模拟星上复杂的振动环境,也不能充分地检验振动补偿系统的性能。本文采用国外实测卫星平台的振动数据模拟星上随机振动,并根据随机振动的特点定义了补偿效率因子ηR以描述补偿效果。搭建平台并开展了验证实验,取得了良好的效果。

1 卫星平台随机振动仿真

最早的星上在轨振动测试在1984年由美国的NASA对Landsat-4测试完成[5],这是人们首次通过实验测得卫星在轨时的功率谱密度(Power Spectral Density,PSD)曲线。1990年,ESA对半导体激光星间通信实验卫星(Emiconductor Intersatellite Laser Experiment,SILEX)进行了在轨微振动测试实验[6],并模拟了一条可以模拟该卫星平台微振动的功率谱密度曲线。

上述两条曲线有很多共同点,如微振动的振幅随着频率的增加而减小,这使得对光通信影响较大的振动主要集中在低频。

有研究者根据上述两次实验测得的功率谱,设计了相应的滤波器,采用对高斯白噪声进行滤波的方法,将频域的功率谱密度转化为时域信号,并以此作为模拟卫星平台微振动的控制信号[7]。该方法得到的信号可以更加真实、有效地模拟星上随机振动,本文中将采用该方法模拟卫星平台上的随机振动。SILEX平台的模拟振动信号如图1(a)所示,其功率谱估计与ESA实测功率谱的对比如图1(b)所示,证明该信号能够模拟SILEX平台的振动。图1(a)中看似连续的振动信号实际由多个离散的点连接而成,点间的时间间隔为0.10s,这组信号可以模拟频率在10 Hz以内的随机振动。实验中制作了多组不同频率的振动信号以研究补偿系统的性能。

2 随机振动补偿原理

补偿系统的控制信号流图如图2所示,图中GDriver(s)、GPZT(s)、GAngle(s)、HFeed Back(s)等参数代表和快速偏转镜相关的传递函数,GDetect(s)是探测环节的传递函数,GComp(s)是补偿环节传递函数,U(t)、α(t)分别是系统的输入与输出,θ(t)是外界扰动,即卫星平台的微振动。

图中与快速偏转镜闭环控制环节的传递函数:

整个系统的传递关系:

本实验中为简化起见,将输入U置零,则在理想情况下,振动完全得到补偿,则输出α也应为零,则补偿环节的传递函数:

其中快速偏转镜的传递函数可以通过实验测量[8,9],CCD探测振动的传递函数则需通过实验标定测得。

3 模拟实验结果及分析

在搭建模拟实验平台的过程中,选用快速偏转镜(Fast Steering Mirror,FSM)作为模拟卫星平台振动、补偿振动的器件,该器件采用压电陶瓷驱动,精度可达µrad量级,响应频率接近1 k Hz[10],被广泛应用于自适应光学等领域[11]。目前,面阵CCD已被用于动态目标跟踪[12,13],本实验平台采用高帧频面阵CCD探测振动,其最高帧频可达200 Hz以上,能够探测低频振动。

振动补偿模拟实验平台的原理框图如图3所示,实验装置实物图如图4所示,由计算机1控制快速偏转镜FSM1模拟卫星平台微振动引起的光束抖动,CCD通过探测光斑的位置移动,实时将角度数据传输至计算机2计算振动引起的光束偏角,并控制快速偏转镜FSM2补偿振动引起的光束抖动。

以往的振动补偿仿真实验中,定义振动补偿因子η评价补偿效果[3]:

振动补偿因子可以有效地评价补偿系统对正弦振动的补偿效果,但是对于随机振动补偿效率的评价则显得不够客观。考虑到随机振动的特点,补偿效率评价应从统计的角度出发,并能够体现补偿系统所抵消的振动能量,本文提出采用随机振动补偿效率因子ηR作为评价手段。

式中:iC、iA分别为补偿前后的光斑偏移量。

当时,补偿效率ηR=100%,振动完全被补偿,光斑位置保持在平衡位置固定不动。这是一种理想情况,对于随机振动,不可能做到完全补偿;

当时,补偿效率ηR=%0,表明补偿系统对所探测到的振动信号无动于衷,振动没有得到任何补偿;

当时,补偿效率ηR<%0,表明补偿系统虽然探测到了振动信号,但是做出了不恰当的补偿,不但没有补偿原有振动的影响,反而推波助澜,恶化了不良影响。这是在构建补偿系统时最不希望看到的结果。

随机振动补偿实验的结果如图5所示,FSM1模拟SILEX平台上的随机振动,将CCD的采样频率设置为100 Hz,理论上该系统能补偿50 Hz以内的振动。图中可见,随着振动频率的增加,补偿效率呈下降趋势。当CCD采样频率高于振动频率三倍以上时,补偿效率高于50%,效果较好。而当CCD采样频率仅为振动频率二倍至三倍时,补偿效果不够理想。可见,提高探测器件的采样频率,可以有效地提高补偿效率,增大补偿带宽。

4 结论

卫星平台的随机振动是影响星间链路通信质量的重要因素,有必要构建振动补偿系统。补偿系统补偿的对象是卫星平台振动中的低频分量。针对随机振动补偿的特殊性,应采用补偿效率因子ηR描述补偿效果。实验结果表明,当探测器的采样频率高于振动频率三倍时,补偿系统能够取得较好的效果。下一步可通过提高探测器件采样频率来提高补偿效率。

摘要：卫星平台随机振动对星间激光链路通信质量产生增大误码率等不良影响。为了有效地补偿随机振动,提高通信质量,本文设计并搭建了一套用于模拟卫星平台振动与补偿的实验系统。该系统通过面阵CCD探测光束偏转角度,采用压电陶瓷致动的快速偏转镜模拟、补偿卫星平台的随机振动。针对随机振动补偿实验的特殊性,文中定义了补偿效率因子ηR以描述补偿效果。实验结果表明,该系统对卫星平台振动中的低频分量具有一定的补偿能力,当CCD采样频率高于振动频率三倍以上时,该系统对低频随机振动的补偿效果显著。

随机振动 第2篇

随机振动对星载TDICCD影响分析

研究了星载TDICCD的.成像机理,着重分析卫星扰动引起TDICCD随机振动,造成图像质量下降,给出了典型的振动传递率和可行的振动抑制方法.为了提高TDICCD减振的能力,对TDICCD实施隔振和阻尼抑振控制技术.同时分析了振动控制的基本原理和设计原则,给出了相关的设计.

作 者：陈丁跃 周仁魁 李英才  作者单位：陈丁跃(长安大学,西安,710064;中国科学院西安光机所,西安,710068)

周仁魁,李英才(中国科学院西安光机所,西安,710068)

刊 名：光子学报  ISTIC PKU英文刊名：ACTA PHOTONICA SINICA 年，卷(期)： 33(10) 分类号：V445.8 关键词：随机振动   TDICCD相机   卫星   振动抑制  

高精度星敏感器随机振动分析 第3篇

星敏感器是卫星和航天器姿态控制系统中不可或缺的组成部分.它利用对星点目标的探测、与星表的匹配计算确定卫星和飞行器的姿态, 是航天设备中非常重要的组成部分.航天器发射过程是一个非常严酷的环境, 这就要求像星敏感器这样的高精度仪器必须要有足够的稳定性.在航天器发射过程中, 星敏感器受发动机工作噪声及气动力激振等因素影响, 将产生随机振动载荷, 这种振动载荷极有可能对星敏感器产生有害影响, 振动过程中的最大应力可能使星敏感器上的关键器件瞬时失灵, 或者导致连接螺栓的松脱甚至断裂.因而需要通过随机振动分析与试验考核星敏感器承受随机振动环境的能力[1], 提高产品的可靠性.

1 振动理论

振动可分为确定性振动和非确定性振动两类.诸如简谐振动、周期振动和非周期振动这些可以用时间的确定函数来表达, 都属于确定性振动.而另一类不能用确定函数来描述, 瞬时值不可知的振动属于非确定性振动.这里主要研究非确定性振动即随机振动.随机振动由于不能用确定性函数描述运动规律, 只能通过概率、统计方法表述随机过程的重要特征.这种振动不可预测, 在相同的条件下也不重复, 具有明确的随机性.随机振动分析是一种采用功率谱密度作为输入的谱分析, 是一种确定响应出现特定值的概率大小的分析方法.

测得的星敏感器随机振动信号往往是加速度—时间信号, 一般不可能通过求得随机信号的傅氏变换而了解随机过程的频率组成.解决这个问题的方法是利用相关函数和功率谱密度[2].

相关分析是在时域中研究随机信号的统计特性的一种方法, 相关函数有自相关函数和互相关函数, 这里以自相关函数为例.一个时间函数x (t) 的自相关函数Rx (τ) 定义为

$\begin{array}{l}R{}_x(\tau )=E[x(t)x(t+\tau )]=\underset{t\to \alpha }{{\displaystyle \lim }}\frac{1}{{\rm T}}{\displaystyle {\int }_{-{\rm T}/2}^{{\rm T}/2}}\\ x(t)x(t+\tau )\text{d}t(1)\end{array}$

式中, E[x (t) x (t+τ) ]表示求x (t) 和x (t+τ) 的数学期望.

相关分析在频域的对应分析方法是功率谱密度分析.自相关函数的傅氏变换是自功率谱密度函数.假设x (f) 是x (t) 的傅氏变换, G (f) 是单边功率谱密度, S (f) 是双边功率谱密度.则有

S (f) =2∫∞0R (t) cos2πfτdτ (2)

$G(f)=\{\begin{array}{l}S(f),f=0\\ 2S(f),f>0\end{array}(3)$

x (t) 的均方值$\overline{x{}^2}={\displaystyle {\int }_{-\infty }^{\infty }S}(f)\text{d}f$, 可见$\overline{x{}^2}$的值可由功率谱密度曲线S (f) ～f下的面积给出.对自相关函数做傅氏分析不会丢失过程的频域信息.对于一个特定的时间域信号可以通过式 (1) ～式 (3) 得到信号的频域信息.

2 星敏感器随机振动分析

随机振动分析可以分为2种方法, 一种方法是如文献[1]中给出的利用经典公式进行估算;另一种方法是利用有限元软件对随机振动试验进行仿真模拟, 这种方法运算量大, 但精度较之前一种高很多, 这也是现在常用的一种方法.星敏感器属于复杂结构, 分析结果的精确度由结构模型的简化程度、网格划分的精度和加载条件等决定.下面主要介绍利用Patran和Nastran有限元软件对星敏感器进行随机振动分析.

2.1 建立有限元分析模型

文中用于分析的星敏感器三维图如图1所示[3,4,5].对星敏感器结构划分有限元网格, 在建立有限元模型时由于结构本身比较复杂, 在不影响分析结果的前提下做了适当的简化, 比如透镜和镜框全接触和透镜曲面以平面代替等.建立的有限元模型如图2所示, 共包含15 791个单元和18 849个节点.

2.2 功率谱密度的确定

功率谱密度实际上就是将对时间域的振动描述转化为对频率域的振动描述.转化过程中利用式 (1) ～式 (3) 可以得出功率谱密度的表达式.利用公式的计算过程比较复杂, 在实际操作中可以简化:将所测得的加速度—时间信号分割成许多小的频率段, 每段只含有一条很窄的频率段的加速度信息, 当每个信号进入频率段后被平方, 然后与前面的信号平均得到这个频率段的平均值, 再除以频率段的宽度就得到了功率谱密度[6].该分析中得到功率谱密度分布见表1所示.

为了计算频率为20 Hz和2 000 Hz时的功率谱密度, 假设频率和功率谱密度有如下对应关系如表2所示, 则满足式 (4)

A=10log (p2/p1) /log2 (f2/f1) (4)

代入数据可以得到20 Hz时对应的功率谱密度为0.004 243 g2/Hz, 20 000 Hz时对应的功率谱密度为0.001 6 g2/Hz.功率谱密度图纵轴单位为g2/Hz, 横轴单位为g2/Hz, g2是频率带加速度均方值.因此计算功率谱密度曲线下的面积所得的值代表了整个频率范围内的加速度均方值, 与第一节中的定义吻合, 将该值开根号所得到平均值即为均方根加速度值[6].利用这个方法计算得到的本次分析均方根加速度值为6.59 g.

2.3 随机振动分析

随机振动用于考虑结构在某种统计规律分布的载荷作用下的随机响应.随机振动分析的基本步骤可分为2步, 第1步是进行频率响应分析;第2步是用MSC.Random 进行随机响应分析.这里对星敏感器沿X轴和Y轴2个方向施加载荷, 分别得出水平和竖直放置时的随机振动响应.

高精度星敏感器随机振动被认为是线性、平稳、各态历经的随机物理过程.对星敏感器随机振动结果分析后提取了有限元模型中的3个特殊点:反应输入加速度载荷的安装环上的固定节点24 677;遮光罩上加速度响应RMS值最大点即离固定点最远的位于Y向边缘的节点15 593;镜框上加速度响应RMS值最大点即离固定点最远的位于Y向边缘的节点4 073 .

加速度响应如图3所示, 其中图3a是点24 677的加速度响应曲线;图3b是X向加载时节点15 593和4 073的加速度响应曲线;图3c是Y向加载时节点15 593和4 073的加速度响应曲线.从图3中可以看出输入加速度与给定分析条件是吻合的, 加速度响应RMS值最大点都位于节点4 073处.加速度放大倍数在X向加载时为2.227, Y向加载时为1.023, 都在高精度星敏感器要求的安全裕度内.

星敏感器应力主要集中在光学镜头靠近固定点镜筒直径骤降的地方和安装环连接安装平台的固定点, 因此应力分布结果给出了这部分单元在频率变化范围内的应力变化曲线, 如图4所示, 图4a为X向加载时应力集中部部分单元应力分布曲线, 图4b为Y向加载时应力集中部部分单元应力分布曲线.从结果可以知道应力最大值在安装环连接安装平台的固定点处.从图4中可以看出, X向和Y向加载两种情况下的应力值都是非常小的, 完全能满足设计要求.

3 结 束 语

从振动理论入手, 阐述了高精度星敏感器随机振动的相关分析过程, 包括功率谱密度的确定和振动分析的结果评价.由分析结果可知, 航天器发射过程中产生的随机载荷对某高精度星敏感器的影响非常小.不论是X向安装还是Y向安装, 星敏感器的最大加速度响应放大倍数都在安全裕度内, 应力集中处应力值都远小于材料的许用值.接下来的工作是结合星敏感器其他静态和动态指标, 在满足要求的同时优化星敏感器的结构.

参考文献

[1] 余成武, 卢欣. 有限元分析在航天器产品设计中的应用[J]. 空间控制技术与应用, 2008, 34 (4) :28-32.

[2] 李德葆, 陆秋海. 工程振动试验分析[M].北京:清华大学出版社, 186-193.

[3] Jeffrey W Percival, Kenneth H Nordsieck, Kurt P Jaehnig. The ST5000: a high-precision star tracker and attitude determination system[C] //Proceedings of SPIE, 2008, 7010: 70104H-1-70104H-6.

[4] DAN Michaels, JAMES Speed. New Ball Aerospace star tracker achieves high tracking accuracy for a moving star field[C]//Proceedings of SPIE, 2004, 5430: 43-52.

[5] ROELOF W H, VAN Bezooijen. SIRTF autonomous star tracker[C]//Proceedings of SPIE, 2003, 4850: 108-121.

[6] 李蓓蓓. 振动分析的有效工具——功率谱密度[J]. 包装工程, 2004, 25 (3) :46-50.

参考文献

[1]余成武, 卢欣.有限元分析在航天器产品设计中的应用[J].空间控制技术与应用, 2008, 34 (4) :28-32.

[2]李德葆, 陆秋海.工程振动试验分析[M].北京:清华大学出版社, 186-193.

[3]Jeffrey W Percival, Kenneth H Nordsieck, Kurt PJaehnig.The ST5000:a high-precision star tracker andattitude determination system[C]//Proceedings ofSPIE, 2008, 7010:70104H-1-70104H-6.

[4]DAN Michaels, JAMES Speed.New Ball Aerospace startracker achieves high tracking accuracy for a moving starfield[C]//Proceedings of SPIE, 2004, 5430:43-52.

[5]ROELOF W H, VAN Bezooijen.SIRTF autonomousstar tracker[C]//Proceedings of SPIE, 2003, 4850:108-121.

随机振动 第4篇

随机振动载荷动力学等效的一种工程实现方法

利用实验室振动加载技术来等效模拟实际工况的随机振动载荷,以研究结构的振动环境适应性,是工程界通用的做法.文中讨论了随机振动载荷动力学等效的过程与方法,针对小阻尼稀疏模态结构,给出了基于结构振动响应等效的随机振动载荷等效关系的工程应用表达式和评价等效载荷对疲劳损伤影响的`方法.简支梁数值模拟表明,通过文中方法设计的等效随机振动载荷,不仅可以获得变化不大的结构响应,而且疲劳损伤等效结果也可以满足工程应用精度要求.

作 者：朱学旺 刘青林 ZHU Xue-wang LIU Qing-lin 作者单位：中国工程物理研究院总体工程研究所,四川绵阳,621900刊 名：实验力学 ISTIC PKU英文刊名：JOURNAL OF EXPERIMENTAL MECHANICS年，卷(期)：22(6)分类号：O324 O346.2关键词：随机振动 振动载荷 动力学等效 振动响应 小阻尼结构 稀疏模态

随机振动 第5篇

导引头工作的力学环境比较复杂, 为了检验其在真实环境下的工作情况, 保证其对环境的适应性和使用的可靠性, 国家军用标准中规定军品必须进行一定的力学环境试验, 随机振动试验是比较重要的一种。因为振动是产品失效的主要环境因素之一, 大多数振动环境是随机振动的[1]。随机振动试验可以有效地暴露产品的早期故障, 提高军品的使用可靠性。在产品的设计和改进阶段一般需要对其进行随机振动分析, 这样可以及早的找出产品的重要缺陷, 对结构设计进行优化, 避免生产浪费, 缩短研制周期, 降低成本。

随机振动分析也称功率谱密度分析 (PSD) , 它属于一种定性分析。功率谱密度是结构对随机动力载荷响应的概率统计, 其原始数学模型是以概率理论为基础的, 与其它分析不同, 在力学上不是一个能够定量分析的问题, 但即使这样, 还是能够从PSD分析中获得一些定性的数据, 如1σ或者3σ位移、速度、加速度以及单元的应力结果, 这里的1σ和3σ响应值就是概率统计中正态分布下的均方根响应值小于该值的出现概率分别为68.27%和99.74%。

在产品的随机振动试验中, 由于受振动试验设备条件的限制, 一般是三个轴分别加载, 但是产品在使用环境中可能三个轴方向上同时受到振动, 因此, 用ANSYS分析产品在三个轴向同时加载的情况下的响应有很大的意义。

1 建立有限元模型

尽管ANSYS的建模技术日益强大, 但是和专业的三维建模软件Pro E相比, 其效率还是相差很多。并且几乎所有的几何模型都是在Pro E中绘制的, 因此, 直接把Pro E中的几何模型导入ANSYS中将大大提高建模效率。

1.1 导入几何模型

把Pro E中的几何模型导入ANSYS中需要对几何模型做一些简化和修改, 主要是去掉对模型分析结果影响不大的细节, 例如小螺钉、垫片等小零件以及螺纹孔、倒角等特征。

1.2 定义几何模型的有限元元素

由于导入的模型是装配体, 因此相互接触的零件之间要定义接触类型, 如果没有特殊要求, 可以接受ANSYS默认定义的接触类型。

对每个零件分配预先定义好的材料, 本文模型中主要有两种材料铝合金和合金钢。

2 模态分析

模态分析是动力学分析的基础, 模态分析用来确定结构的振动特性。由于结构的振动特性决定结构对于各种动力载荷的响应情况, 所以在准备进行其它动力分析之前首先要进行模态分析。

模态分析的理论和方法在这里就不再赘述。边界条件设为位标器下底面为固定约束。模态分析提取了前20阶的固有频率和振型, 本文只给出模型前6阶的固有频率及其前4阶的振型。

3 随机振动分析

在模态分析的基础上, 就可以对产品进行随机振动分析了。

3.1 加载功率谱密度函数

随机振动分

析很重要的步骤就是功率谱密度函数的加载。设计任务书给出的随机振动条件如图2所示, 图谱所示为宽带随机振动加窄带随机振动组合的类型。

在ANSYS里输入随机振动条件时, 需要知道每个频率拐点的谱值。因此还需计算频率在120Hz和180Hz时的谱值。图谱所示曲线坐标为对数坐标, 各频率点的谱值计算方法见式 (1) [3]。

式中, Di为第i点频率的谱值 (g2/Hz) ;m为每倍频程斜率 (-6d B/Oct) ;fi为第i点频率值 (Hz) 。

根据式 (1) , 可以计算出频率在120Hz和180Hz时的谱值分别为0.025g2/Hz和0.011g2/Hz。

3.2 计算分析结果

首先, 把功率谱密度函数分别单独加载到模型的X、Y、Z三个轴上, 计算模型的响应。得出模型的应力云图如图3、4、5所示。然后, 把功率谱密度函数同时加载到模型的三个轴上, 得出模型的应力云图如图6所示。

从图中可以看出, 在对X、Y、Z三个轴分别加载时, 模型所对应的最大应力分别为71.3MPa、81.2MPa和90.0MPa。当对模型三个轴同时加载时, 模型最大的应力为138.6MPa。

可见, 对模型三个轴同时加载时比单个轴独立加载时所受的应力要大。因此, 在设计阶段对模型作随机振动分析时, 最好用三个轴同时加载的情况。本例中模型三个轴同时加载的最大应力值远小于材料的抗拉强度540MPa。设计结果满足振动条件要求。

4 结束语

本文以导引头位标器为研究对象, 利用有限元分析软件AN-SYS, 根据导引头工作的力学环境条件, 对导引头模型进行了模态分析和随机振动分析。对单轴独立加载和三轴同时加载的随机振动分析进行了对比。结果表明, 用三个轴同时加载的情况更符合产品实际使用情况。此方法也说明ANSYS可以仿真出由于受振动试验设备条件的限制而不能作的试验, 对产品的设计和改进有一定的参考意义。

参考文献

[1]王敏, 王小静.弹载天线随机振动的仿真分析[J].制导与引信.2007;28 (2) :57-60.

[2]刘衍平, 高新霞, 张刘斗等.电子设备机柜随机振动试验的数值模拟[J].塑性工程学报.2007;14 (4) :151-155.

[3]崔冠宇, 李亮.电子设备减振设计的随机振动分析.遥测遥控[J].2011;32 (3) :68-73.

[4]范文杰.星载电子设备宽频随机振动响应分析[J].电子机械工程.2010;26 (4) :5-17.

随机振动 第6篇

Ye Mao,Xiao Feng,Mai Zhendong,et al.Evolutionary random vibration response of a coupled vehicle--bridge system with sprung masses.Mechanics in Engineering,2015,37(4):475-480

引言

随着经济的快速发展，高架桥在现代城市交通中扮演着越来越重要的角色.桥梁结构越来越柔长，而车辆却向高速化、重载化方向发展，这使得车桥耦合振动越来越明显.因此，车桥耦合系统的减振日益成为人们关注的重点.文献[1,2]采用双臂弯矩控制系统的主动控制方法，研究移动载荷作用下桥梁的主动最优控制.文献[3,4,5,6]采用在桥梁上悬挂TMD(turred mass damper)的半主动控制方法研究车桥耦合系统的减振控制.不论采用哪种减震形式，在设计计算时，大多忽略桥面不平度的影响,或者由给定的轨道不平顺功率谱,通过三角级数叠加法[3]模拟桥梁不平顺.大多数研究者在建立车-桥--TMD或弹簧-质量系统体系的运动方程时，单独对TMD或弹簧-质量系统建立振动方程，将悬挂点的相容条件加入到耦合方程中.

文献[7]通过在梁上加弹簧-质量系统进行减振，并获得良好的效果.本文从演变随机过程的一般理论出发，以桥面不平度功率谱密度函数作为系统的能量输入.在建立车-桥-弹簧-质量系统运动微分方程时，不单独对弹簧-质量系统建立运动方程，而是将弹簧-质量系统的运动考虑到桥梁的模态中去.提取桥梁-弹簧-质量系统的模态，将该模态代入车桥耦合系统，计算系统的随机响应.

1 振动方程的建立

1.1 车--桥耦合振动方程

带弹簧-质量系统的车桥耦合系统模型如图1所示.简支梁在x1,x2，···,xn处分别布置n个弹簧-质量系统，其质量分别为mT1,mT2,mTn；刚度分别为kT1,kT2,kTn；桥梁总长L，桥梁单位长度质量ρA，桥梁抗弯刚度EI，单位长度阻尼系数c.车辆简化为两自由度模型，如图1所示.u(x,t),u1(t),u2(t)分别表示桥梁、车辆两个质量m1和m2的竖向位移响应.车辆在桥梁上的位置为ξ(t).η(ξ)为桥面不平度值.不单独对弹簧-质量系统建立运动方程，在车体质量与桥梁质量相比很小的情况下，根据达朗贝尔原理，车-桥耦合系统的动力方程为

根据模态分析法，梁的变形u(x,t)可以写成如下形式

式中，ϕi(x)是连续梁桥的第i阶模态，qi(t)是对应于第i阶模态的正则坐标.

将式(2)代入式(1)，利用模态函数关于质量和刚度的正交性，得如下微分方程组

式中mb=ρAL.

方程组(3)在时域内有n+2个变量，令q=[q1q2···qnu1u2]T，可将方程组(3)写成矩阵形式

式中，M,C,K分别为耦合系统的质量、阻尼和刚度阵；f1(t)是车辆自重产生的确定性激励，f2(t)是桥面不平顺产生的随机激励，其表达式如下

其中b=[0,0,···,0,k2]T.

1.2 桥--弹簧--质量系统的运动方程与模态求解

对于具有n个弹簧-质量系统的桥梁，kT为悬挂弹簧的刚度，悬挂质量为mT,梁截面高度为d,弹性模量为E，梁单位长度质量为ρA，截面惯性矩为I，梁的总质量为mb,mTi,¨zTi,kTi,zTi分别表示第i个弹簧-质量系统的质量，加速度，刚度，位移.YTi表示悬挂弹簧-质量振子处梁的位移.建立运动方程为

弹簧振子处的变形协调条件为

其中，分别表示梁的位移、转角、曲率.表示弹簧振子i的左侧梁截面，表示弹簧振子i的右侧梁截面.

桥梁与质量之间的平衡条件为

简支梁的边界条件为

根据协调条件，平衡条件和边界条件求解式(2)即可得到桥梁与弹簧-质量系统耦合的模态.求解的过程文献[8]已详细介绍.将求得的模态带入到车-桥耦合系统中即可求解振动方程.

为验证本文程序的正确性，采用文献[8]中的仅带一个悬挂弹簧-质量的两端简支的梁的模型参数,全长L=1 m，弹簧-质量系统安装在离左侧3/4处，kT=1.904 28×105EI/L3,mT=0.2ρAL=3.077 5 kg,E=2.069×1011N/m2,d=0.05 m,ρA=15.387 5 kg/m,I=3.067 96×10-7m4,mb=ρAL=15.387 5 kg得出振型结果如图2所示.图2与文献[8]所得完全一致.由此可以验证本文采用的求解模态方法与程序的正确性.

2 随机响应

方程(4)的载荷项由两部分组成，即确定性激励与随机激励，二者之间相互独立，可分别计算.耦合系统的总响应为两者叠加.本文仅考虑由桥面不平顺引发的随机响应.即，令f1(t)=0.

令，则方程(4)可写为

式中

式(9)计算由桥面不平顺引发的耦合系统响应桥面不平度η(ξ)是随自变量ξ变化的零均值均匀随机场，功率谱密度为S(ω).

零初始条件下，式(9)的解可表示为

式中，Φ(t,ε)为式(9)的齐次方程的转换矩阵.

令η(t)=η[ξ(t)]，η(t)的相关函数可表示为

式中，ξ(t)表示车辆在桥上的位置.当车辆匀速行驶时，令ξ=vt，代入式(11)可得

可得到B(t)的相关矩阵

状态空间向量κ的相关矩阵为

由式(12)可知，式(13)可以改写为

式中

为H(ω,t)的共轭转置矩阵.

对照式(10)，从式(16)可看出，当ω为定值时，H(ω,t)就是零初始条件下，原系统在确定性激励作用下的响应.可以通过数值方法求解，如Runge--Kutta法.

由式(13)的结果，可得桥梁截面垂向位移均方值响应，即

设Hi(ω,t)表示矩阵H(ω,t)第i行，式(15)可改写为

根据式(18)，式(17)变为如下形式

由此，得到桥梁截面位移的演变功率谱密度函数为

同样可得车体质点的位移响应均方值

相应的演变功率谱密度函数为

桥面不平度的谱密度函数一般由位移谱给出，其速度谱与位移谱有下列换算关系

因此

车体及桥梁截面处的垂向加速度均方响应如下

响应的加速度演变功率谱密度函数分别为

3 数值算例

选取文献[9]桥梁参数，具体为：桥梁总长L=25 m，弹性模量E=2.87 GPa，截面惯性矩I=2.9 m[4]，桥梁单位长度质量ρA=2 303 kg/m，阻尼系数取c=4 681 N·s/m.车辆参数取文献[5]，具体为：车体质量m1=38 500 kg，悬挂体系质量m2=8 660 kg，其刚度与阻尼分别为k1=3.82×10[5]N/m,c1=5.07×10[6]N·s/m,k2=1.96×10[5]N/m.车辆行驶速度V=40 m/s.

没有悬挂弹簧-质量系统的“空”梁的前4阶频率为[30.020 1,120.080 6,270.181 3,480.322 2]rad/s.在桥梁中点悬挂一个弹簧-质量系统.其质量和刚度与桥梁参数的关系为

算例中取µ=0.008，ν=3.

此时，整个梁与弹簧-质量系统体系的前4阶频率为[24.618 2,70.981 6,120.080 6,272.903 9]rad/s.其前4阶模态如图3所示.

桥梁不平度η(ξ)为零均值均匀随机场，根据《机械振动——道路路面谱测量数据报告》标准(GB/T 7031-2005)和文献[10]，采用下列桥面不平度功率谱密度函数在时域内的表达式

式中,V表示车辆行驶速度，ω为圆频率(rad/s)，ω0为参考频率，ω0=2πn0V,n0为参考空间频率，n0=0.1m-1；ω的取值范围为[2πnminV,2πnmaxV]，其中nmin=0.011,nmax=2.83；本文采用A级桥面等级，则桥面平整度系数为Sr(n0)=16×10-6m-3.

在有弹簧-质量系统与无弹簧-质量系统两种工况下，桥梁中点的随机响应分别如图4所示.从图4中可以明显看出：梁上悬挂弹簧-质量之后，桥梁中点的位移与加速度随机响应都有大幅下降.

质量比µ是弹簧-质量系统的一个重要参数,工程实践中一般取[0.005,0.02][11].为研究质量比µ对桥梁中点随机响应的影响，分别取不同质量比，其他参数同算例.在不同质量比下，桥梁中点随机响应的峰值如图5所示.

从图5(a)中可以看出：桥梁中点的位移随机响应随质量比的逐步增大有一个先下降后上升的趋势在质量比µ=0.016时，桥梁中点的随机响应最小图5(b)中可以看出随质量比的增大，桥梁中点的加速度响应逐渐减小.

在工程应用中也有在桥梁上布置多个弹簧-质量系统的情形.根据前面的讨论,此时取弹簧-质量系统的参数为：µ=0.016，ν=3.

图6给出了3种工况下桥梁中点的随机响应时程曲线.工况(1)：单个弹簧-质量系统，安装位置在跨中；工况(2)：两个弹簧-质量系统，安装位置分别在距左端1/4和3/4处.工况(3):3个弹簧-质量系统，安装位置分别在距左端1/4,2/4和3/4处

从图6可以看出：3种工况下，工况(2)时，桥梁中点的位移和加速度均方根响应峰值最小，而工况(3)时其位移峰值响应最大，工况(1)时其加速度最大.由此可见，在桥梁上安装弹簧-质量系统时，对桥梁随机响应的控制，并不是弹簧-质量系统个数越多就减振效果越好.因此，建议在实际工程中适当位置布置弹簧-质量系统，使桥梁的随机响应降低

4 结论

将弹簧-质量系统的振动考虑到桥梁的模态中去，提取桥梁-弹簧-质量系统的模态，以桥面不平度功率谱为输入，采用Runge--Kutta法计算车桥耦合系统的演变随机响应.最后，得出弹簧-质量系统可以减少桥梁随机响应，适当的弹簧-质量系统质量比参数、安装位置可以大幅降低桥梁的随机响应

参考文献

[1]宋娟,熊冶平,张玉双.移动荷载下桥梁的随机性主动最优控制.东南大学学报,2001,31(5A):52-54

[2] 宋娟,熊冶平,王锡平.移动荷载下桥梁的振动主动控制.振动工程学报,2001,14(S):83-85

[3] 袁小钦,刘习军,张素侠.MR-TMD减振系统对连续梁桥振动控制研究.振动与冲击,2012,31(20):153-157

[4] 肖新标,沈火明.移动荷载作用下的桥梁振动及TMD控制.振动与冲击,2005,24(2):58-61

[5] 黎明安,崔凯,雷霜.车桥耦合振动系统的半主动控制法.应用力学学报,2011,28(4):376-380

[6] 郭文华,路萍.TMD对高速列车通过简支箱梁时的振动控制研究.振动与冲击,2008,27(11):42-45

[7] 周翔,李俊峰.利用弹簧一质量系统抑制梁上任意点的振动.力学与实践,2007,29(2):18-20

[8] Wu JS,Chou HM.A new approach for determining the natural frequencies and mode shapes of a uniform beam carrying any number of sprung masses.Journal of Sound and Vibration,1999,220(3):340-357

[9] 殷新锋,蔡春声,刘扬.基于面接触的车桥耦合振动研究.振动工程学报,2012,25(3):244-252

[10] 方同.工程随机振动.北京:国防工业出版社,1995

随机振动 第7篇

20世纪80年代以来,农用车逐渐增多,而农村路面状况未得到改善,驾驶员的乘坐舒适性亟待提高。因此,各种线性弹性元件与特定机构配合实现非线性弹性特性的座椅悬架相继问世,极大地改善了农用车驾驶员的乘坐条件。非线性座椅悬架的出现解决了线性悬架的大部分问题。目前,国内外学者对非线性座椅进行的研究,普遍是通过物理样机试验进行分析[1,2,3],或通过软件建模进行仿真分析[4]。通过建立能够反映客观实际的动力学模型,利用Matlab软件计算出农用车受路面随机激励时人体的振动响应量,对座椅的减振特性进行研究,既不受实验条件的限制,又可缩短研究周期,并且提高计算精度,为进一步改善座椅的乘坐舒适性提供了理论基础。

1 人-椅-路面动力学数学模型

1.1 数学模型的建立

本文以X型非线性座椅悬架为原型,它由拉、压缩弹簧、空气阻尼器和平行四杆机构组成,同时采用预紧结构,如图1所示。左右各一对铰接的刚性杆作为支撑体,每对杆的一头与座椅铰接,另一头可沿座椅上的限位滑槽滑动,其与拉伸弹簧、压缩弹簧和空气阻尼器共同起到减缓座椅振动的作用

为了便于研究,根据文献[5],将整个人—椅—路面振动系统简化为6个自由度的振动模型。人体分为头部、胸部器官、骨盆和腿部4个部分,农用车及座椅分别简化为1个自由度,其振动系统如图2所示。

系统的运动微分方程组为

式中m1,m2,m3,m4,m5,m6—农用车、座椅、腿部、骨

盆、胸部器官和头部质量(kg);

c1,c2,c3,c4,c5,c6—农用车、座椅、腿部、骨盆、

胸部器官和头部阻尼系数(N·s/m);

k1,k2,k3,k4,k5,k6—农用车、座椅、腿部、骨盆、

胸部器官和头部刚度系数(N/m);z0—路面的位移(m);

z1,z2—农用车和座椅相对于地面的垂向位移(m);

z3,z4,z5,z6—人体各部分相对于地面的垂向位移(m)。

其中,座椅的刚度为

式中l—座椅支撑杆长度(m);

θ0—拉杆的初始角(rad);

θ—拉杆的运动角(rad);

α—长度选择系数;

k2—主压缩弹簧刚度(N/m);

k″2—拉伸弹簧刚度(N/m)。

根据座椅的结构,将式(2)中θ的转换成y的函数。y是座椅相对于农用车的位移,y=z2-z1,将其用泰勒级数展开,取其前3项,则座椅的刚度可表示为

农用车行驶的多为农村道路,地面不平度系数不能参照国家8级路面不平度系数。根据英国、德国、法国等国家的测量和统计资料,将农村农业地面不平度的自功率谱密度函数用以下形式的指数函数近似描述[6]即

在进行随机振动分析时,根据此函数不同形式的输入得出不同的输出。如在求解时,考虑座椅悬架的非线性项,则式(1)属于非线性随机振动方程。但是,目前尚无普遍有效的求解方法,许多问题只能在特殊情况下求解,大量问题仍停留在近似解决阶段。目前,常用的几种求解非线性系统受随机激励时响应的数值方法有FPK法、摄动法和统计线性化方法。

根据非线性随机振动理论,对采用FPK法、统计线性化方法与将系统线性化所计算出来的均方响应进行分析比较,结果表明:在计算均方值时考虑座椅非线性项的一、二、三项与只考虑一项时的结果误差是1.362 8%,此结果为将系统线性化提供了依据。为了简化计算,本文将图2所示的非线性振动系统线性化,再运用多自由度线性系统的随机振动求解方法对式(1)进行求解。

1.2 数学模型的验证

根据式(1)建立的数学模型采用多自由度系统随机振动的求解方法,利用Matlab软件对求解过程进行编程。将计算出来的理论传递率曲线与物理样机实验的实测传递率曲线进行对比,如图3所示。

从图3可以看出:两条曲线基本吻合,这说明利用Matlab软件来求解式(1)所建立的数学模型能够比较真实地反映实际情况;同时也说明将非线性系统线性化是可行的,可以根据此模型对座椅的各个参数进行优化。

2 参数优化

根据座椅的减振原理和相关理论分析表明,主压缩弹簧的刚度k2、拉伸弹簧(两根)的刚度k″2和阻尼器的阻尼系数C直接决定了座椅的减振性能,所以选定k2、k″2和C为研究因素。通过Matlab软件的计算技术,采用三因子二次回归正交设计进行参数优化,对数学模型进行初步的计算后,确定各因素的水平,如表所示

根据表1,需对式(1)所建立的数学模型进行15次计算试验,将计算出来的数据以ISO 2631振动规范曲线为评价标准。本文采取其中的疲劳—降低功效界限作为评价座椅振动舒适性的标准。

驾驶员主要承受垂向振动,而其身体部分对频率为4～8Hz范围内垂向振动最敏感。对比15组计算出来的1～80Hz频率范围内人体各部位加速度均方值,结果表明:当主压缩弹簧的刚度k2=29 823N/m、拉伸弹簧的刚度k″2=8 511N/m、阻尼系数C=1 082N·s/m时,在4～8Hz频率范围内的加速度值比其他参数组合所计算出来的加速度值小,而且4～8Hz频率范围内头部及胸部器官的加速度值基本位于振动规范8h曲线以下,系统固有频率处的加速度值也都在1h曲线以下,如图4所示。根据ISO 2631振动标准:若在不属于人体敏感频带范围内有一部分振动加速度值超出允许界限,而低于另一稍放宽的允许界限,则仍可认为符合标准。因此,认为采用上述参数的座椅悬架隔振效果最好。根据强度理论对按上述参数设计的弹簧进行强度校核,结果表明符合强度要求。

3 对比试验

为验证优化参数后座椅悬架的减振性能,进行了农村路面行驶对比试验。试验在某农场进行,试验路面起伏较大,相当于较差的农村土路,试验区段长为1 000m。选择两种典型的农用车及拖拉机对X型悬架座椅进行试验:一是东方红902拖拉机;二是菲亚特拖拉机。试验结果表明,将优化参数后的X型悬架座椅依次装于以上两种车型,并以同样的速度沿着同样的路线试验3次,并用DH 5938振动测试分析系统测量人体各部位垂直方向加速度加权均方根值,结果显示:在人体敏感频带4～8Hz范围内,人体各部位的加速度均方值均位于振动规范8h曲线以下。从上述试验结果得出,优化参数后X型悬架座椅可以达到较好的减振效果,是一种适合农用车辆的减振座椅。

4 结论

1)Matlab软件计算随机振动下人体各部位的加速度响应与物理样机实测数据结果相吻合,表明把随机振动下人—椅—路面组成的非线性系统等价为线性系统正确,可用于对座椅的各参数进行优化,缩短研究周期,提高计算精度。

2)通过对X型座椅悬架的分析可知:当座椅的主压缩弹簧刚度k2=29 8231N/m、拉伸弹簧的刚度k″2=8 511N/m、阻尼系数C=1 082N·s/m时,减振效果好,人体在敏感频率4～8Hz范围内的加速度值最小座椅的乘坐舒适性最佳

3)采用优化参数后,X型座椅悬架的农用车辆人体乘坐舒适性明显得到提高,减振能力平均提高40%左右

摘要：针对农村路面崎岖不平、农用车驾驶员对乘坐舒适性要求提高等问题,以减振理论为基础,通过非线性随机振动理论,建立农用车多自由度“人—椅—路面”非线性振动系统数学模型,利用Matlab软件进行人体加速度响应的分析求解,并通过ISO2631国际标准的振动舒适性曲线进行评价,采用正交回归设计方法优化农用车X型非线性座椅的各参数,为改进座椅的乘坐舒适性提供理论基础。路面试验表明:采用优化参数后,农用车X型悬架座椅的人体乘坐舒适性明显得到提高,减振能力平均提高40%左右。

关键词：农用车,座椅,随机振动,乘坐舒适性,参数优化

参考文献

[1]李素云,杨坚.手扶拖拉机X型座椅垂直方向振动的虚拟优化研究[J].拖拉机与农用运输车,2005,2(1):29-31.

[2]方瑞华.汽车空气悬架非线性振动理论和试验[J].农业机械学报,2007,38(7):13-15.

[3]张润生,牛毅.东方红902履带拖拉机乘坐舒适性研究[J].拖拉机与农用运输车,2000(4):12-14.

[4]舒明霏,陈思忠.基于M atlab的越野汽车非线性悬架系统仿真[J].计算机仿真,2006,23(12):216-219.

[5]Cho-Chung Liang,Chi-Feng Chiang.A study on biody-namic models of seated hum an subjects exposed to verticalvibration[J].International Journal of Industrial Ergonom ic,2006,36:869-890.

随机振动 第8篇

作为数学分析的一个重要分支,分数阶微积分的历史可以追溯到300年前。近几十年来,它被广泛应用于反常扩散、信号处理与控制、流体力学、图像处理、黏弹性阻尼器及分形理论等的研究中。特别地,在工程实际中,分数阶导数可以利用较少的参数来模拟黏弹性材料依赖于频率的结构阻尼[1,2,3],使得它被引入到结构动力学的研究中。其中,以分数阶导数来描述黏弹性结构阻尼的分数阶微分方程是这一研究的重点。

随机微分方程是在确定性结构方面,Laplace变换[1,2]和Fourier变换[4]被用于线性分数阶微分方程的求解。而对于非线性情形,至今没有非常有效的求解办法。文献[5]中,平均法被应用于非线性分数阶微分方程,分析了振动结构的振幅响应。实际中,往往会有随机荷载作用于受黏弹性阻尼的振动结构,形成含有分数阶导数阻尼的随机振动结构及描述相应运动的分数阶随机微分方程。文献[6]利用随机平均法来得到响应的近似稳态概率密度。分数阶微分方程的理论和研究方法可见文献[7]。

相对于解析方法,数值方法是分数阶微分方程研究的另一个重点和难点。分数阶导数所具有的历史依赖性及全域相关性导致了它在数值计算方面的复杂性。在数值模拟中,随着时间历程的增加,计算量呈现指数增长。此外,分数阶导数的定义还不完善,至今没有一个统一的定义形式被广泛接受,这也增加了分数阶导数数值计算方法开发的难度。现阶段,分数阶导数的数值算法主要包括:有限差分法、级数逼近法、有限元法、无网格方法等[8,9,10]。随机模拟本身要求进行大量的样本计算,所以对于分数阶随机微分方程的研究要求一些高效的数值算法被提出。

本文旨在给出针对黏弹性阻尼下随机振动结构的数值模拟方法。对于以分数阶导数描述的结构阻尼,利用经典的Riemann-Liouvill定义与Grunwell-Letnikov定义,给出对历史数据进行合理截断的理论分析,在一定意义上减低分数阶导数对历史的依赖性,“缩短”它的记忆性,从而提高数值模拟的效率。对于所提出的数值方法,给出数值算例,通过模拟结果与解析结果的对比来验证其有效性。

1 分数阶导数的定义

分数阶微积分自诞生以来就有多种定义形式,其中分数阶导数普遍采用的主要有三种定义:Riemann-Liouvill(RL)定义、Caputo(C)定义和Grunwald-Letnikov(GL)定义。这三种定义各有定义式,并分别有其适用性,在某些前提下,定义之间是等价的。

1.1三种定义的表达式

1.1.1 RL定义

$D{}^{\alpha }x(t)=\frac{1}{\Gamma (m-\alpha )}\frac{\text{d}{}^m}{\text{d}t{}^m}{\displaystyle {\int }_0^t(}t-\tau ){}^{m-\alpha -1}x(\tau )\text{d}\tau$; m-1≤α<m (1)

式(1)中m是不小于α的最小整数。RL定义可以理解为先对函数x(t)进行m-α次积分,然后再进行m次微分。这个定义具有良好的分析性质,比较适用于对随机过程的分析。

1.1.2 C定义

$D{}^{\alpha }x(t)=\frac{1}{\Gamma (m-\alpha )}\frac{\text{d}{}^m}{\text{d}t{}^m}{\displaystyle {\int }_0^t(}t-\tau ){}^{m-\alpha -1}x{}^m(\tau )\text{d}\tau$;

m-1≤α<m (2)

从式(2)式中可以发现,C定义先对函数x(t)进行m次微分,再进行m-α次积分。C定义下,各初始状态的物理意义与一般整数阶微分下相同,使得该定义可以很好的与实际工程应用背景相结合。

1.1.3 GL定义

$\begin{array}{l}D{}^{\alpha }x(t)={\displaystyle \sum _{k=0}^{m-1}\frac{x{}^k(0)t{}^{-\alpha +k}}{\Gamma (-\alpha +k+1)}}+\\ \frac{1}{\Gamma (m-\alpha )}{\displaystyle {\int }_0^t(}t-\tau ){}^{m-\alpha -1}x{}^m(\tau )\text{d}\tau ;\end{array}$

m-1≤α<m (3)

GL定义可以理解为RL定义的离散形式,可以用来对RL定义进行数值计算。

1.2三种定义的比较

RL定义及C定义都可以看成是对GL定义的一种改进,在阶数α为负实数和正整数时,这两种定义是等价的。而GL定义作为一种级数定义,可以用来对分数阶导数进行数值计算。在函数x(t)具有m阶连续导数的情况下,RL定义与GL定义是等价的,若无上述条件,RL定义是对GL定义的扩充。从定义式也可以明显看出,当初始条件满足xk(0)=0; k=0,1,2,…,m-1的条件下,C定义与GL定义等价。当阶数α是整数时,各种定义均退化为一般意义下的整数阶导数。

2 一种分数阶导数阻尼下随机振动结构的

数值模拟方法

考虑如下的单自由度随机振动结构

$\begin{array}{l}\ddot{X}(t)+\epsilon cD{}^{\alpha }X(t)+g(X,\dot{X})=\\ f(t)+\epsilon {}^{1/2}{\displaystyle \sum _{k=1}^{{\rm M}}h}{}_k(X,\dot{X})W{}_k(t)(4)\end{array}$

X(t)是随机振动振子的位移;DαX(t)表示X(t)对时间t的α阶导数,它描述了作用于振子的黏弹性结构阻尼。当阶数α等于1时,它表示振子所受的线性阻尼,当阶数α等于0时,这一项表示振动结构的刚度。实际中,阶数α是一个介于0到1之间或者1到2之间的任意实数,分别对应欠阻尼和过阻尼的情形。g(X,$\dot{X}$)是X和$\dot{X}$的函数,表示了振子受到的黏弹性阻尼之外的阻尼及回复力。f(t)是确定性的外界激励。hk(X,$\dot{X}$;k=1,2,…,M是(X,$\dot{X}$)的函数,表征了乘性或加性噪声Wk(t); k=1,2,…,M的作用幅度,其中随机荷载Wk(t)用高斯白噪声来模拟,其统计特性为

$\begin{array}{l}E[W{}_k(t)]=0(5)\\ E[W{}_k(t)W{}_1(t+\tau )]=2D{}_{kl}\delta (\tau )(6)\end{array}$

这里Kkl表示白噪声的强度。

式(4)中的分数阶导数DαX(t)采用RL定义

$\begin{array}{l}D{}^{\alpha }X(t)=\\ \{\begin{array}{l}\frac{1}{\Gamma (1-\alpha )}\frac{\text{d}}{\text{d}t}{\displaystyle {\int }_0^t(}t-\tau ){}^{-\alpha }X(\tau )\text{d}\tau ;0<\alpha <1\\ \frac{1}{\Gamma (2-\alpha )}\frac{\text{d}{}^2}{\text{d}t{}^2}{\displaystyle {\int }_0^t(}t-\tau ){}^{1-\alpha }X(\tau )\text{d}\tau ;0<\alpha <2\end{array}(7)\end{array}$

如前所述,该定义便于进行理论分析,但是并不便于进行数值计算,所以首先对RL定义的分数阶导数利用GL定义进行处理。振子的位移X(t)在弱阻尼的情形下可以被认为具有2阶以上连续导数,则式(7)对应的GL定义为

$\begin{array}{l}D{}^{\alpha }X(t)=\\ \{\begin{array}{l}\frac{X(0)t{}^{-\alpha }}{\Gamma (1-\alpha )}+\frac{1}{\Gamma (1-\alpha )}{\displaystyle {\int }_0^t(}t-\tau ){}^{-\alpha }X{}^{(1)}(\tau )\text{d}\tau ;\\ 0<\alpha <1\\ \frac{X(0)t{}^{-\alpha }}{\Gamma (1-\alpha )}+\frac{X{}^{(1)}(0)t{}^{-\alpha +1}}{\Gamma (2-\alpha )}\frac{1}{\Gamma (2-\alpha )}\times \\ {\displaystyle {\int }_0^t(}t-\tau ){}^{1-\alpha }X{}^{(2)}(\tau )\text{d}\tau ;0<\alpha <2\end{array}(8)\end{array}$

进行分部积分,得到

$D{}^{\alpha }X(t)=\frac{1}{\Gamma (-\alpha )}{\displaystyle {\int }_0^{t}X}(t-\tau )(\tau ){}^{-\alpha -1}\text{d}\tau (9)$

对其中的积分利用积分定义进行离散

$\begin{array}{l}D{}^{\alpha }X(t)=\frac{1}{\Gamma (-\alpha )}\underset{h\to 0}{{\displaystyle \lim }}{\displaystyle \sum _{k=0}^{{\rm N}}X}(t-kh)(kh){}^{-\alpha -1}h=\\ \underset{h\to 0}{{\displaystyle \lim }}h{}^{-\alpha }{\displaystyle \sum _{k=0}^{{\rm N}}X}(t-kh)\frac{k{}^{-\alpha -1}}{\Gamma (-\alpha )}(10)\end{array}$

由函数的定义

$\begin{array}{l}\Gamma (-\alpha )=\\ \underset{k\to \infty }{{\displaystyle \lim }}\frac{k!k{}^{-\alpha }}{(-\alpha )(-\alpha +1(-\alpha +2)\cdots (-\alpha +k)}(11)\end{array}$

可以得到

$\begin{array}{l}\frac{k{}^{-\alpha -1}}{\Gamma (-\alpha )}\to (-1){}^k\frac{\alpha (\alpha -1)(\alpha -2)\cdots (\alpha -k+1)}{k!}\stackrel{\Delta }{{\displaystyle =}}\\ (-1){}^k\left(\begin{array}{l}\alpha \\ k\end{array}\right);k\to \infty (12)\end{array}$

则有

$\begin{array}{l}D{}^{\alpha }X(t)=\underset{h\to 0}{{\displaystyle \lim }}h{}^{-\alpha }{\displaystyle \sum _{k=0}^{{\rm N}}(}-1){}^k\left(\begin{array}{l}\alpha \\ k\end{array}\right)X(t-kh)=\\ \underset{h\to 0}{{\displaystyle \lim }}h{}^{-\alpha }{\displaystyle \sum _{k=0}^{{\rm N}}G}L{}_{k}X(t-kh)(13)\end{array}$

这里

$GL{}_k=(-1){}^k\left(\begin{array}{l}\alpha \\ k\end{array}\right);k=0,1,2,\cdots ,{\rm N}(14)$

称为GL系数。从定义式可以得到,由于

$\left(\begin{array}{l}\text{α}\\ 0\end{array}\right)=1$

,则GL0=1;又根据

$(-1){}^k\left(\begin{array}{l}\alpha \\ k\end{array}\right)=\frac{\Gamma (k-\alpha )}{\Gamma (-\alpha )\Gamma (k+1)}(15)$

并根据Γ函数的递推关系

$\Gamma (z)=(z-1)\Gamma (z-1)(16)$

得到

$\begin{array}{l}GL{}_k=\frac{\Gamma (k-\alpha )}{\Gamma (-\alpha )\Gamma (k+1)}=\frac{(k-\alpha -1)}{k}\times \\ \frac{\Gamma (k-\alpha -1)}{\Gamma (-\alpha )\Gamma (k)}=\frac{k-\alpha -1}{k}GLk-1(17)\end{array}$

分数阶导数与整数阶导数最显著的区别在于它们在计算时是否需要考虑已经过去的状态。整数阶导数的计算只涉及到与所考虑时刻相关的一个小临域内的函数取值,与其他时刻没有关系,即所谓的“无记忆性”或“局部性”;而分数阶导数依赖于整个过去状态,具有“记忆性”或“整体性”,在计算时必须考虑到全部历史数据。显然,在GL系数的递推关系式(17)中,$\frac{k-\alpha -1}{k}<1$,所以GL系数逐项递减,说明了分数阶导数在某一时刻的值对于历史数据的依赖性随着时间的远离而减弱,即某时刻的分数阶导数值对相对较近的历史时刻的依赖性较大,而对于相对较早的历史时刻的依赖性较小。所以在实际的计算中,为了提高计算的效率,可以对分数阶导数对历史时刻的记忆性进行合理的截断,只考虑相对较近的时刻,忽略更早时刻的影响,利用较短的记忆来估算某一时刻的分数阶导数值。

对于给定的阶数α,可以计算相应的GL系数序列GLk。例如对于α=0.8,可以得到如图1所示的GL系数序列。从图中可以看出,当k>5时,GL系数已经达到10-2量级,而当k>20时,GL系数可以达到10-4量级,所以可以认为更早时刻的函数取值对此时刻分数阶导数的计算没有影响,在计算时忽略不计。

基于以上的截断思想,对于式(13)可以得到以下的分数阶导数的近似计算公式

$\begin{array}{l}D{}^{\alpha }X(t)=\underset{h\to 0}{{\displaystyle \lim }}h{}^{-\alpha }{\displaystyle \sum _{k=0}^{{\rm N}{}_{tru}}G}L{}_{k}X(t-kh)\approx \\ \underset{h\to 0}{{\displaystyle \lim }}h{}^{-\alpha }{\displaystyle \sum _{k=0}^{{\rm N}{}_{tru}}G}L{}_{k}X(t-kh)(18)\end{array}$

这里用Ntru来表示截断后所剩的项数,Ntru具体的取值取决于阶数α。至此,分数阶导数的近似数值计算已经解决,代回到式(4),可以得到分数阶随机微分方程的数值计算方法。这里可以采用的迭代方法很多,本文采用随机Runge-Kutta方法,对式(4)进行Monte Carlo模拟。

3 数值算例

为了验证所提出数值方法的有效性,将其应用于受分数阶导数阻尼的线性随机振动结构

$X(t)+cD{}^{\alpha }X(t)+\omega {}_0{}^{2}X(t)=W(t)(19)$

式(19)中c为阻尼系数,ω0为振子的自然频率,W(t)是噪声强度为D的高斯白噪声。为了得到响应的近似解析解,首先引入变换

$\{\begin{array}{l}X(t)=A(t)\cos \Phi (t)\\ \dot{X}(t)=-A(t)\Omega (A,\Phi )\sin \Phi (t)\end{array}(20)$

将对原响应(X,$\dot{X}$)的讨论转换为对慢变变量(A,θ)的讨论。将变换式(20)代入式(19),得到

$\{\begin{array}{l}\frac{\text{d}A}{\text{d}t}=\frac{cA\Omega \sin \Phi }{g(A)}D{}^{\alpha }(A\text{c}\text{o}\text{s}\Phi )-\frac{A\Omega \sin \Phi }{g(A)}W(t)\\ \frac{\text{d}\theta }{\text{d}t}=\frac{cA\Omega \cos \Phi }{g(A)}D{}^{\alpha }(A\text{c}\text{o}\text{s}\Phi )-\frac{\Omega \cos \Phi }{g(A)}W(t)\end{array}(21)$

式(21)中回复力g(X)=ω02X。运用标准的随机平均法,得到漂移系数m(A)和扩散系数σ(A)。

$\{\begin{array}{l}m(A)=c\left(\frac{A\Omega \sin \Phi D{}^{\alpha }(A\text{c}\text{o}\text{s}\Phi )}{g(A)}\right){}_{\Phi }+\\ D\left(\frac{\partial G{}_1}{\partial A}G{}_1+\frac{\partial G{}_2}{\partial \theta }G{}_2\right){}_{\Phi }\\ \sigma {}^2(A)=2D(G{}_1{}^2){}_{\Phi }\end{array}(22)$

式(22)中$G{}_1=-\frac{A\Omega \sin \Phi }{g(A)}‚G{}_2=-\frac{\Omega \cos \Phi }{g(A)}$。为了得到漂移系数m(A),第一项的确定性平均涉及到分数阶导数的积分,这里需要应用Fresnel积分的渐近积分公式

$\begin{array}{l}{\displaystyle {\int }_0^t\frac{\cos (\overline{\Omega }\tau )}{\tau {}^{\alpha }}}\text{d}\tau =\\ \overline{\Omega }{}^{\alpha -1}\left[\Gamma (1-\alpha )\sin \frac{\alpha \text{π}}{2}+\frac{\sin (\overline{\Omega }t)}{(\overline{\Omega }t){}^{\alpha }}+0[(\overline{\Omega }t){}^{-\alpha -1}]\right](23)\\ {\displaystyle {\int }_0^t\frac{\sin (\overline{\Omega }\tau )}{\tau {}^{\alpha }}}\text{d}\tau =\\ \overline{\Omega }{}^{\alpha -1}\left[\Gamma (1-\alpha )\cos \frac{\alpha \text{π}}{2}-\frac{\cos (\overline{\Omega }t)}{(\overline{\Omega }t){}^{\alpha }}+0[(\overline{\Omega }t){}^{-\alpha -1}]\right](24)\end{array}$

经过计算,即得到关于振幅A的Ito随机微分方程

$\text{d}A=m(A)\text{d}t+\sigma (A)\text{d}B(t)(25)$

及描述其概率密度演化的FPK方程

$\frac{\partial p}{\partial t}=-\frac{\partial }{\partial A}[m(A)p]+\frac{1}{2}\frac{\partial {}^2}{\partial A{}^2}[\sigma {}^2(A)p](26)$

它的稳态解为

$p(A)=\frac{C}{\sigma {}^2(A)}\exp \left\{{\displaystyle {\int }_0^A\frac{2m(u)}{\sigma {}^2(u)}}\text{d}u\right\}(27)$

对于振动结构式(19),运用上述的标准随机平均法,得到随机振动振子振幅响应的稳态概率密度

$p{}_{st}(A)=\frac{C\omega {}_0{}^{2}A}{D}\exp \left\{-\frac{c\omega {}^{1+\alpha }\sin (\alpha \text{π}/2)}{2D}A{}^2\right\}(28)$

即振幅A服从Rayleigh分布,其中C为归一化常数,根据式(28)得到振子的位移和速度的联合稳态概率密度

$\begin{array}{l}p{}_{st}(X,\dot{X})=\\ \frac{C\omega {}_0}{2\pi D}\exp \left\{-\frac{c\omega {}_0{}^{1+\alpha }\sin (\alpha \pi /2)}{2D}(X{}^2+\dot{X}{}^2/\omega {}_0{}^2)\right\}(29)\end{array}$

相应可以得到关于位移和速度的稳态概率密度

$p{}_{st}(X)=\frac{C\omega {}_0}{2\pi D}\exp \left\{-\frac{c\omega {}_0{}^{1+\alpha }\sin (\alpha \text{π}/2)}{2D}X{}^2\right\}(30)$

$p{}_{st}(\dot{X})=\frac{C\omega {}_0}{2\pi D}\exp \left\{-\frac{c\omega {}_0{}^{1+\alpha }\sin (\alpha \text{π}/2)}{2D\omega {}_0{}^2}\dot{X}{}^2\right\}(31)$

即振子的位移和速度分别服从正态分布${\rm N}\left(0‚\frac{D}{c\omega {}_0{}^{1+\alpha }\sin (\alpha \pi /2)}\right)$和${\rm N}\left(0‚\frac{D\omega {}_0{}^2}{c\omega {}_0{}^{1+\alpha }\sin (\alpha \pi /2)}\right)$。图2中,在一组给定参数下,当分数阶阶数取0.5时,用近似解析结果对所提出的随机模拟方法进行验证。图2(a)给出了此时的GL系数,从第40项开始,GL系数减小到10-2量级,可以将这之前的历史数据舍弃,即Ntru取50,只考虑它之后的项对某一时刻分数阶导数的影响。图2(b),(c)及(d)给出了振子稳态概率密度的近似解析解和数值解的比较。其中实线为式(30)及式(31)给出的解析结果,“o”表示数值模拟的结果。从图中可以看出,对于典型的欠阻尼随机振动结构,提出的数值方法可以很好的对其进行数值模拟,而且,利用GL系数对分数阶导数的记忆性进行有效截断之后,可以大大的提高随机模拟的效率。图3中,分数阶导数的阶数取1.5,类似图3的分析,得到此时的近似解析解与数值结果的对比,两者很好的吻合,说明在过阻尼情形,所给出的数值计算方法也是有效的。

为了进一步验证数值算法的有效性,图4给出了不同的分数阶阶数时随机振动位移响应的方差,其中参数为D=0.01,c=0.05,ω0=2,实线表示解析结果,“o”是对应阶数下的模拟结果,可以看出,数值解与近似解析解可以很好的吻合,说明所给出的随机模拟方法是有效的。

4 结论

提出了一种针对受分数阶阻尼的随机振动结构的数值计算方法。在对比了分数阶导数各种定义的基础上,使用在理论中具有明显优势的RL定义来模拟随机振动结构中的分数阶阻尼,并用GL定义来对其进行离散处理,以便进行数值计算。给出了GL系数的递推公式,并指出了合理的阶段阶数,使得截断后的分数阶导数值既能较好的反应分数阶导数对于历史数据的记忆性,又能提高数值模拟的效率。数值算例说明所提出的数值方法能够很好地吻合随机振动结构响应的解析解,包括振幅、位移和速度的稳态概率密度及位移的方差,从而说明了该数值方法的有效性。

参考文献

[1]Bagley R L,Torvik P J.Fractional calculus-a different approach tothe analysis of viscoelastically damped structures.AIAA Journal,1983;21(5):741-748

[2]Bagley R L,Torvik P J.Fractional calculus in the transient analysisof viscoelastically damped structures.AIAA Journal,1985;23(6):741-748

[3]Koeller R C.Application of fractional calculus to the theory of vis-coelasticity.ASME Journal of Applied Mechanics,1984;51(2):299-307

[4]Gaul L,Klein P,Kemple S.Impulse response function of an oscilla-tor with fractional derivative in damping description.Mechanics Re-search Communications,1989;16(5):297-305

[5]Wahi P,Chatterjee A.Averaging for oscillations with light fractionalorder damping.Proceedings of ASME 2003Design Engineering Tech-nical Conferences and Computers and Information in Engineering Con-ference;Chicago,2003:721-727

[6]Huang Z L,Jin X L.Response and stability of a SDOF strongly non-linear stochastic system with light damping modeled by a fractional de-rivative.Journal of Sound and Vibration,2009;319(3-5):1121-1135

[7] Podlubny I. Fractional Differential Equations. London: Academic Press 1999

[8]Diethelm K,Ford N J,Freed A D et al.Algorithms for the fractionalcalculus:A selection of numerical methods.Computer Methods inApplied Mechanics and Engineering,2005;194(6-8):743-773

[9]Diethelm K,Ford N J,Freed A D.A Predictor-Corrector Approachfor the Numerical Solution of Fractional Differential Equations.Non-linear Dynamics,2002;29(1-4):3-22