分数阶黏弹性振子系统响应的代理模型分析

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

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

随着国家经济的不断进步，科学技术的不断发展，工程车辆已经普遍应用于工程应用中，为工程的实施带来了便利。但是对于工程车辆来说，由于其本身工作环境的特殊性，使得其在特殊工况下工作时会产生剧烈的振动及噪声，这不仅会对车辆的使用寿命有所损耗，还会对驾驶员的身心健康有所影响。研究表明，黏弹性材料作为阻尼所构建的黏弹性缓冲结构由于其廉价的制作成本和高效的减振缓冲能力，使得其被广泛地应用于车辆的减振控制中，以达到更好的缓冲减振效果[1]。同时由于分数阶建模优于整数阶建模的非局部性，使得其仅需较少的本构参数就能更加准确的描述黏弹性材料的遗传特性，能够更好地反映黏弹性材料的力学特性。所以，分数阶黏弹性建模被广泛应用于工程车辆的研究中[2-3]。

随着分数阶微分建模的不断发展，分数阶微积分方程的解法也层出不穷。目前已有一些较为经典的解析解算法，如Laplace变换、Fourier变换等[4]。但是由于分数阶黏弹性系统计算的复杂性，使得其解析解难以求得，因此对于分数阶微分方程数值解的研究已成为学者们研究的热点，常见的数值算法有：有限差分法[5]、谱方法[6]、正交多项式[7]逼近算法、Adomian分解法[8]和小波分析法[9]等。Ezz-Eldien[10]首次尝试用谱方法求解FLPMs，结果表明：该方法能够用较小的移位多项式获得较高的精度。Mashali-Firouzi等[11]提出了一种用于求解非线性多项分数阶初值问题的多步自适应伪谱方法，该方法求解精度高且适合大范围的计算。Kim等[12]运用移位雅可比频谱伽辽金法来求解分数阶微分方程的初值问题，并验证了多阶线性和非线性分数阶微积分方程解的收敛速度，但是求解过程较为繁琐，在实际的工程问题中，需要耗费较高的成本。

近年来，数学模型与计算机模拟的结合已成为解决物理问题的方法之一[13]。但是由于计算机强大的模拟能力使得其运算成本较高，进而增大了实际工程中的成本预算。代理模型的出现很好地解决了这个问题，用代理模型或者元模型替代昂贵的工程问题，从而降低工程中的运算成本。代理模型是工程优化中常见的一种方法和手段，其思想是通过一定的方式和手段来处理优化过程中复杂的计算问题。代理模型具有较好的仿真模拟能力，被广泛应用于各种工程问题中，包括优化设计、可靠性分析、基于可靠性分析的设计优化等[14]。

鉴于此，以标准性固体本构的分数阶黏弹性振子(fractional viscoelastic oscillators of solid state constitutive, SFVEO)系统作为研究对象，首先给出一种基于分数阶导数之间关系的数值差分格式，并分析在单位冲激激励下系统参数对位移响应的影响。最后，为了解决分数阶黏弹性系统优化计算的复杂性问题，用代理模型对系统的位移响应进行代理。研究成果可为分数阶黏弹性振子系统响应分析及工程设计提供理论参考。

1 SFVEO模型建立及求解

1.1 SFVEO模型的建立

如图1所示，SFVEO系统由两个理想胡可弹性元k1、k2及分数阶黏弹性元组成，其中，c为黏弹系数，α、β为分数阶阶次。分数阶黏弹性元与弹簧k1为串联结构，与弹簧k2为并联结构。

σ、ε分别为黏弹性材料的应力和应变

黏弹性材料的分数阶本构关系可表示为

Pσ=Qε

(1)

式(1)中：P、Q为表达黏弹性材料力学特性的分数阶本构算子，其表达式分别为

(2)

(3)

Caputo定义下的分数阶微分算子，其定义为[15]

(4)

式(4)中：α为分数阶阶次；f(t)为求导函数；t为自变量；n为常数；Γ为Gamma函数；τ为中间变量。

取n=1,p0=1，p1=c/k1+k2，β0=0，β1=β，m=1，α0=0，α1=α，q0=k1k2/k1+k2，q1=k1c/k1+k2，则该黏弹性系统的本构关系为

σ+p1Dβσ=q0ε+q1Dαε

(5)

当系统在外界激励的作用下，系统发生形变，为便于描述该系统的力学特性，现给系统外接一个质量为m的物体，如图2所示，当系统受到的外部激励所施加的力为f时，系统发生形变所产的阻尼力大小为fd。

图2 SFVEO模型受力图

根据经典振动力学方程有

(6)

(7)

由动力学方程得

(8)

式(8)中：l为该分数阶黏弹性系统得长度；x为系统振动时外接物产生的位移。

将式(6)～式(8)代入式(5)中，并结合分数阶微分算子的性质得

(9)

给定初值，将其转化为数学模型为

(10)

式(10)中：a1=p1；a2=1；a3=q1ζ/m；a4=q0ζ/m；b1=1/m；b2=p1/m；c1、c2、c3为常数。

1.2 SFVEO模型的求解

(11)

式(11)中：ωj为算术矩阵，无实际意义；j表示矩阵的行列数；h为步长；n=「α⎤，当初值为0时，分数阶Caputo导数定义与分数阶G-L导数相互等价。

用分数阶G-L算子代替分数阶Caputo算子，则方程的解为

(12)

式(12)中：Pn(x)为外界激励和分数阶导数项的非零初值项；(FM)t为外界激励的零初值项。

2 单位冲激激励下系统阻尼特性分析

当ξ=0.5，ζ=0.05，ωn1=ωn2=10，β=0.1时，分数阶阶次α变化对系统响应曲线的影响如图3所示。

图3 阶次α变化对系统响应曲线的影响

当α取值小于0.5时，如图3(b)所示，随着α取值的不断增大，响应曲线的第一幅值随之减小，响应曲线的衰减速度不断加快；当α取值大于0.5时，如图3(c)所示，随着α取值的不断增大，响应曲线的衰减速度越来越快，但与图3(b)相比，α取值越大，对系统的影响越不明显。由图3(a)可知，对于SFVEO模型，当调整模型参数时，会导致系统的响应曲线出现过阻尼的现象。

当ξ=0.5，ζ=0.05，ωn1=ωn2=10时，阶次β变化对系统响应曲线的影响如图4所示。由图4可知，阶次α和β的取值之间存在相互抑制的关系，当α=0.8时，β=0.5时出现过阻尼现象；当α=0.6时，β在0.4处就已经出现过阻尼现象。即α的取值越大时，β的取值范围就越大。当阶次α的取值固定时，随着β取值的不断增大，系统的振幅衰减速度逐渐减弱，且对系统位移响应影响较小。

图4 阶次β变化对系统响应曲线的影响

当ζ=0.05，ωn1=ωn2=10，α=0.7，β=0.2，阻尼比ξ=0、0.5、1时，系统响应曲线的变化如图5所示。

图5 阻尼比ξ=0、0.5、1时系统响应曲线的变化

随着阻尼比取值的不断增大，系统响应曲线的第一幅值呈下降趋势，衰减速度不断加快。

当ξ=0.5，ωn1=ωn2=10，α=0.7，β=0.2，几何参数ζ=0.05、0.5、5时，系统响应曲线的变化如图6所示。

图6 几何参数ζ=0.05、0.5、5时系统响应曲线的变化

单位冲激激励下，随着几何参数的不断增大，系统响应曲线的幅值较大，衰减速度较快。

当ξ=0.5，ζ=0.05，α=0.7，β=0.2，固有频率ωn1=10、50、100，ωn2=10、50、100时，系统响应曲线的变化如图7所示。

图7 ωn1=10、50、100，ωn2=10、50、100时系统响应曲线的变化

当ωn1处于低频状态时，随着ωn2的不断增大，系统的衰减速度不断增大，在低频段，ωn2对系统响应曲线的幅值及衰减速度的影响较大。随着ωn1取值的不断增大，ωn2对系统响应曲线幅值及衰减速度的影响逐渐减小。当ωn1=10，ωn2处于低频阶段时，响应曲线衰减速度较快，当处于中频或高频阶段时，系统响应曲线的第一幅值和衰减速度较为接近，对系统响应曲线影响较小。当ωn1和ωn2均处于中频或高频阶段时，系统响应曲线衰减不明显。

3 代理模型的构建

3.1 代理模型类型及构建流程

目前，常用的代理模型有Kriging模型、响应面模型、径向基函数模型、正交多项式模型、支持向量回归模型等。Kriging代理模型和响应面代理模型在处理复杂问题时计算时间短且精度高[16]，因此主要对这两个模型进行介绍。

3.1.1 Kriging代理模型

Kriging模型是一种基于统计理论的插值模型，由向量x的多项式模型和向量x的局部偏差组合而成。该模型不仅能够给出未知函数的预测值，还能给出预测值的误差估计，其模型结构为

(13)

3.1.2 多项式响应面模型

多项式响应面模型是一种基于最小二乘法的多项式回归模型，拟合结果能够以显示函数的形式表示出来，便于后续优化设计，且具有较好的连续性和可导性，使用非常便利。在多项式响应面模型中，若x是一个独立因子向量，y为响应值，则存在式(14)所示的关系。

y=f(x)+e

(14)

式(14)中：e为均值为0、标准差为0的正态分布随机误差。

常用的模型有一次多项式响应面模型、二次多项式响应面模型、三次多项式响应面模型等。由于二次多项式响应面模型中无需较多的参数变量且变量间存在交互性，故被广泛应用于求解实际工程问题，其模型表达式为

(15)

代理模型的优化过程大致分为试验设计(design of experiment, DOE)、模型建立及优化分析3部分。其流程图如图8所示。

图8 代理模型构建流程图

由图8可知，模型的建立步骤如下。

步骤1选取所研究系统的参数作为设计变量。

步骤2选取合适的试验设计方法。

步骤3从系统中选取训练样本点和预测样本点。

步骤4选择代理模型。

步骤5通过训练样本的不断训练代理模型，生成训练后的代理模型。

步骤6运用训练后的代理模型预测出预测样本点的预测值。

步骤7对预测样本点的预测值与真实值进行误差分析。

步骤8设置精度要求，若训练后的代理模型达到了精度要求，则输出最终的代理模型；若不符合精度要求，则返回步骤3，扩大样本点的采点数量，从而继续训练代理模型，直至达到精度要求。

3.2 SFVEO系统代理模型的构建

依据代理模型构建步骤，现以SFVEO系统作为研究对象，进行代理模型的构建，构建具体过程如下。

3.2.1 设计变量的选取

设计变量的选取会直接影响优化问题的计算效率及精度，设计变量的数量过多或过少时都会影响代理过程的计算效率及精度。为了能较为准确的模拟出系统的位移响应函数，依据2节中地参数分析进行参数选取，阶次β对系统位移响应反应不敏感，故选取系统频率ωn∈[0,20]、阻尼比ξ∈[0,1]、分数阶阶次α∈[0,1]和几何参数ζ∈[0,1]参数作为设计变量。

3.2.2 试验方法的选取

选取拉丁超立方方法选取样本点，在每个变量的取值范围内选取训练样本点及测试样本点。这里选取50个样本点，且选用Kriging代理模型及二次多项式响应面模型进行代理计算，通过误差分析，从而确定合适的代理模型形式。Kriging模型拟合效果如图9所示。

图9 Kriging模型拟合效果

借助计算机实验的设计与分析(design and analysis of computer experiments，DACE)工具箱构建Kriging模型，设置初样本点个数N为50，选取线性相关模型和二次多项式回归模型，相关函数中的参数θ选取3.152，对SFVEO模型进行代理，当不满足精度要求时选用加点策略，通过增加样本点的个数来提高拟合精度。

由图9可知，随着样本点数量的增加，拟合效果更为准确。同时，选用二次响应面模型对SFVEO进行代理，设置初始样本点个数N为50，选用拉丁超立方设计选取样本点。设置决定系数R2≥0.9时达到精度要求，两种代理模型的误差分析结果如表1所示。

表1 代理模型误差分析对比

由表1可知，在选取较少样本点时，二次响应面模型对SFVEO系统响应函数的代理效果更好，则输出拟合后的SFVEO位移响应函数为

xS=0.009 18-0.000 20ωn1-0.000 14ωn2-

0.001 73ξ-0.001 72α-0.005 93ζ-

0.000 16α2+0.000 04ζ2-0.000 17ωn1ωn2-

0.000 24ωn1ξ-0.000 1ωn1α-0.001 5ωn1ζ-

0.002 42ωn2ξ+0.000 007ωn2α+0.000 01ωn2ζ+

0.000 97ξα+0.001 50ξζ+0.012 30αζ

(16)

综上，相较于Kriging模型而言，二次响应面代理模型在选取少量的样本点时就能对分数阶黏弹性模型进行代理。

4 结论

以SFVEO系统作为研究对象，分析了其在单位冲激激励下参数对系统动态响应的敏感程度，并选用Kriging模型和二次响应面模型对系统的响应函数进行代理，得出如下结论。

(1)α的取值对系统阻尼特性具有较大的影响，但是随着α取值越大，对系统的影响越不明显；β的取值对系统阻尼特性的影响较弱。

(2)随着阻尼比和取值的增大，系统响应曲线的第一幅值呈下降趋势，衰减速度不断加快；随着几何参数ζ取值的增大，响应曲线衰减速度不断加快。

(3)频率ωn1和ωn2在低频阶段对系统阻尼特性的影响较大。当频率ωn1处于高频阶段时，系统出现过阻尼现象。

(4)二次响应面代理模型对分数阶黏弹性系统的代理效果优于Kriging模型。

系统的参数发生变化会导致系统出现过阻尼现象，实际的工程中为了避免这种现象所带来的困扰，在实验设计中需要对系统的刚度-阻尼参数进行优化。用二次响应面代理模型得到的结果可以用于后续的优化计算中，减少优化过程计算复杂度。