基于闭式解算法的粘弹性振子系统阻尼效应分析

孙 宝，张文超，李占龙，秦 园

(1.太原科技大学 应用科学学院， 太原 030024； 2.太原科技大学 机械工程学院，太原 030024)

由于大型工程车辆工作环境的特殊性，使其在作业时会产生严重的振动和噪声，对周边环境造成很多负面的影响。研究发现粘弹性阻尼结构能够很好的解决这一问题，且制造和维修费用较低，所以被广泛的应用到工程车辆的振动控制中。在工程车辆的建模中，主要是通过结构的应力、应变和位移的三变量关系构造系统的基本关系，由于整数阶微分模型的局限性，很多情况下无法准确的表达出粘弹性材料的本构关系及其力学特性，所以分数阶微分已经被广泛应用于粘弹性材料及其缓冲结构动态阻尼分析等方面的建模中[1-2]。

近年来，求解分数阶微分方程的方法层出不穷，但是由于求解分数阶微分方程的解析解难度较大，且一般解的形式都是带有特殊函数的(例如多变量的Mittag-Leffler函数、Green函数等)，这些特殊函数计算起来比较困难，甚至无法得到最终的确切值，所以对分数阶微分方程数值解的研究就尤为重要[3-5]。常见的数值方法有Adomain分解法[6]，变分迭代法[7]，同伦分析法[8]等。此外，任建娅等[9]将Haar小波[10]与算子矩阵思想结合在一起，利用小波矩阵独有的性质求解分数阶微分方程，使得计算更加简便，且精度较高。B.K.Singh等[11]在同伦摄动法的基础上结合Laplace变换提出了一种新的同伦摄动变换法(简称HPTM)，该数值方法所求得的级数形式的解收敛速度快，常用来求解各种工程中的物理模型。Wang Z[12]提出了一种用广义Adams-Bashforth-Moulton方法与线性插值方法相结合的方法来求解一类非线性分数阶微分方程，并验证了其数值格式是有效的。 M.Khan等[13]在同伦分析法的基础上结合Laplace变换思想形成了一种新的同伦分析变换法，该方法在选取合适的参数后，能够用较少的迭代次数求得方程的数值解且精度较高，但是在实际的应用中不易实现。Golbabai A等[14]提出了一种利用Bernoulli多项式的性质求解的方法，但求解过程较为复杂且计算量较大。

本文建立了一种可用于构建大型工程车辆阻尼缓冲结构的Maxwell分数粘弹性振子系统(以下简称MFVEO)。针对MFVEO系统进行建模，根据粘弹性材料的分数阶本构关系及动力学关系，建立分数阶数学模型。采用闭式解算法求解零初值下系统模型的数值解，并对该粘弹性系统的阻尼效应进行分析。

1 MFVEO模型的建立及求解

1.1 预备知识

为了方便对分数阶微分方程的研究，不同的研究者在研究过程中给分数阶微分算子赋予了不同的定义。常见的有Riemann-Liouville、Caputo和Grünwald-Letnikov三种定义，由于Grünwald-Letnikov定义下的微分算子不含积分形式，可写为离散形式，所以在此用到的算子均为Grünwald-Letnikov定义下的分数阶微分算子，定义如下[15]：

定义1：函数f(t)在Grünwald-Letnikov定义下的分数阶微分为

(1)

分数阶微分是对整数阶的一个扩充，所以与整数阶微分算子一样，分数阶微分算子也具有以下性质：

(2)

(3)

Laplace变换是分数阶微分方程中最常用的积分变换，对分数阶常微分方程的求解具有重要的作用，下面对其相关定义进行一些简单的介绍。

定义2：设函数f(t)在[0,+∞]上有定义，对任意一个实数(或复数)s，如果有

(4)

则称F(s)为函数f(t)的Laplace变换，F(s)和f(t)分别为像函数和原函数。

1.2 MFVEO模型的建立

在MFVEO系统中，该阻尼系统是由理想胡克弹性元k和分数阶粘弹性元〈c,α,β〉两部分组合构成的，如图1所示，其中α和β为材料系数，0<α<1,0<β<1。

图1 MFVEO模型结构示意图

设其本构关系为

(5)

给该阻尼系统外接一个质量为m的物体(图2)，当该物体受到一个大小为f的外力时，该阻尼系统发生变形时提供的阻尼力为fd，设S为物体的受力面积，a为该阻尼系统的长度，x为物体受到外力作用时发生的位移，与所受外力方向相反。

图2 MFVEO模型受力分析图

根据物理中的动力学方程及运动学关系，我们可以得到以下关系：

(6)

(7)

(8)

将式(6)、式(7)和式(8)代入式(5)，得：

(9)

根据G-L分数阶微分的性质，上式可写为：

(10)

将式(10)整理后得到MFVEO系统的分数阶微分方程为

(11)

将式(11)进一步简化为

(12)

1.3 MFVEO模型的求解

在零初值的条件下，对方程式(12)的第一个式子两边同时作用于Laplace变换，则得到系统的分数阶传递函数为

(13)

根据式(1)中Grünwald-Letnikov的定义，其离散形式可写为

(14)

其中

(15)

将式(14)代入式(12)的第一个式子中，得：

(16)

将含xt的项提出，整理得：

(17)

所以，由式(17)可知式(12)在零初值条件下的闭式解为

(18)

该方法主要用到了离散的G-L分数阶算子的定义，相对于其他的数值解方法来说，计算量较小，易得出分数阶微分方程的数值解，可用于求解零初值下MFVEO系统的数值解。下面讨论分析不同激励下参数变化对MFVEO阻尼效应的影响。

2 MFVEO系统阻尼效应分析

2.1 正弦激励下系统的阻尼效应分析

现假设系统所受外界激励为正弦激励f=sinωt，其中ω为外力激励频率。取m=80，令系数a1=0.1,a2=0.5,b1=0.012 5,b2=0.001 25,ω=2π，分析MFVEO系统参数改变对系统阻尼效应所带来的影响。

1) 阶数α对系统阻尼效应的影响。固定阶数β的值为0.1且方程系数不变时，分析阶数α取值变化对系统阻尼效应的影响。为了能够更加方便简洁的观察响应曲线的变化，现将其分为α∈(0,0.5)和α∈(0.5,1)两部分来进行讨论，结果如图3所示。

图3 阶数α取不同值时粘弹性振子随时间变化曲线

从图3中可知当β=0.1且方程系数都不发生变化时，当α取值为0.1时，系统响应曲线的波动最为强烈。随着分数阶微分的阶数值α的不断增大，该粘弹性振子振动的幅值不断减小，系统的阻尼效应不断增强。

2) 阶数β对系统阻尼效应的影响。固定阶数α的值为0.3和0.8且方程系数不变时，分析阶数β变化对系统阻尼效应的影响，结果如图4所示。

图4 阶数β取不同值时粘弹性振子随时间变化曲线

由图4可知：无论α取何值，随着β取值的不断增大，该粘弹性振子振动的幅值不断增大，系统的阻尼效应越来越弱。

3) 外力激励频率ω对系统阻尼效应的影响。固定阶数α=0.4，β=0.1，在其他系数不变时，讨论当系统所受外力激励的频率ω取2π，4π和6π时，对系统阻尼效应的影响，结果如图5所示。

图5 频率ω=2π,4π,6π时粘弹性振子随时间变化曲线

由图5可知：随着外力激励频率ω的不断增大，该粘弹性振子的幅值不断减小，系统的阻尼效应不断增大。对比图3和图4，当固定方程的系数不变时，系统的阻尼效应受阶数α的影响较大，受阶数β的影响较小。

2.2 单位冲激激励下系统的阻尼效应分析

当外界激励由正弦激励变为单位冲激激励后，取m=80，令系数为a1=0.1,a2=0.5,b1=0.0125,b2=0.001 25，分析MFVEO系统的阶数发生改变时，对系统阻尼效应所带来的影响。

1) 阶数α对系统阻尼效应的影响。固定阶数β=0.1且方程其他系数不变，分析单位冲激激励下阶数α取值变化对系统阻尼效应的影响。结果如图6所示。

图6 阶数α取不同值时粘弹性振子随时间变化曲线

由图6可知：取β=0.1且方程系数都不发生变化，当阶数α∈(0,0.5)时，随着α取值的不断增大，该粘弹性振子振幅衰减的速度越来越快，系统的阻尼效应不断增强；当阶数α∈(0.5,1)时，随着α取值的不断增大，该粘弹性振子振幅衰减的速度越来越慢，系统的阻尼效应不断减弱。由于粘弹性材料的特性，当α的取值越小，越趋于0时，粘弹性振子的波动曲线越接近于弹簧的振动，即材料的力学特性更趋向于弹性；当α的取值越大，越趋于1时，粘弹性振子的波动曲线的波动变化幅度越小，即材料的力学性能更趋向于粘性。

2) 阶数β对系统阻尼效应的影响。固定阶数α的值为0.3和0.8且方程系数不变时，分析单位冲激激励下，阶数β变化对系统阻尼效应的影响，结果如图7所示。

由图7可知：固定阶数α及方程系数不变时，随着阶数β取值的不断增大，该粘弹性振子振幅衰减的速度越来越快，系统的阻尼效应不断增强。

图7 阶数β取不同值时粘弹性振子随时间变化曲线

3 结论

1) 在正弦函数激励下，MFVEO系统的阻尼效应与阶数α的变化与呈正相关，与阶数β的变化呈负相关；且该系统的阻尼效应对阶数α的变化较敏感，对阶数β的变化不敏感；随着外力频率ω的不断增大，系统的阻尼效应不断增强。

2) 在单位冲激激励下，当阶数α∈(0,0.5)时，随着α取值的不断增大，系统的阻尼效应不断增强；当阶数α∈(0.5,1) 时，随着α取值的不断增大，系统的阻尼效应不断减弱。随着阶数β取值的不断增大，系统的阻尼效应不断增强。

3) 该闭式解算法可用于求解零初值下的分数阶微分方程的数值解。而在实际工程中，系统的初值不一定都为0，所以如何求解初值不为0的分数阶微分方程的解，需进一步探讨和研究。