橡胶衬套高阶分数导数动力学模型的研究*

2019-09-04冯金芝郑松林

汽车工程 2019年8期

高琦，冯金芝，郑松林，林阳

（1.上海理工大学机械工程学院，上海 200093； 2.上海汇众汽车制造有限公司，上海 200122）

前言

橡胶衬套的动态性能影响悬架的动力学特性，它具有较强的非线性黏弹力学特性，受加载幅值、加载频率和工作周期等因素的影响［1］。衬套模型的精度是悬架及整车动力学仿真精度的关键影响因素之一，目前整车及子系统模型动态仿真精度不足很大原因是衬套模型精度不够造成的［2］。因此，建立准确描述衬套动态特性的动力学模型，对提高整车仿真精度、准确预测悬架动态K＆C特性对整车操稳性和平顺性的影响具有重要意义。

目前为止，国内外研究人员提出了多种橡胶元件的动力学模型，为克服Kelvin-Voigt模型、三参数Maxwell模型和由多个Maxwell单元并联得到的广义橡胶衬套动力学模型的不足，近年来提出的分数导数模型因能反映加载历程对橡胶材料动态特性的影响而日益受到关注［3-7］，这类模型均为最高阶次不大于1的低阶分数导数模型，研究人员发现，最高阶次大于1的高阶分数导数模型更接近于描述橡胶材料的物理本质，对黏弹性力学特性的描述更加精确［8-9］，但相关研究文献相对较少，因此，本文中尝试用高阶分数导数模型描述橡胶衬套的黏弹性力学特性，建立橡胶衬套黏弹性高阶分数导数动力学模型。

以某车用悬架摆臂橡胶衬套为研究载体，进行动态加载试验，分析其轴向频率相关性和振幅相关性；建立基于弹性单元、摩擦单元和高阶分数导数黏弹性单元并联的橡胶衬套力学模型，提出一种基于动态邻居广义学习策略粒子群优化（dynamic neighbor and generalized learning particle swarm optimization，DNGL-PSO）算法的高阶分数导数黏弹性单元参数识别方法，并结合试验结果对模型参数进行识别；进一步，为提高模型整体的预测精度和适用性，将黏弹性单元的分数导数系数拟合为加载振幅的函数。通过与试验结果对比，验证所建立的模型能够精确地描述悬架衬套轴向的动态特性。

1 橡胶衬套动态加载试验

基于MTS831衬套试验台进行了某车用悬架摆臂橡胶衬套轴向动态加载试验，激励为振幅和频率一定的正弦信号，试验装置与夹具的安装方式如图1所示。

图1 橡胶衬套轴向动态加载试验

试验对象为轴对称圆柱形橡胶衬套，考虑该衬套在工作过程变形受到约束，这里仅通过动态加载试验研究其在线弹性变形范围内的动态特性。通过分析该悬架衬套的工作特性，进行了频率范围为1～40 Hz、振幅为0.1～2 mm的动态加载试验，得到的橡胶衬套的动刚度和阻尼系数随振幅和加载频率的变化规律如图2和图3所示。

图2 不同加载幅值下衬套动刚度随频率的变化

图3 不同加载幅值下衬套阻尼系数随频率的变化

由图2和图3可以看出：1～40 Hz频率范围内，当振幅一定时，橡胶衬套的动态刚度和阻尼系数都随着频率的增加而增大；随着振幅的增加，橡胶衬套的动刚度显著下降，阻尼系数随着振幅的增加先增大至峰值再呈现出减小的趋势。

2 橡胶衬套动力学建模

所建立的衬套动力学模型由弹性单元、摩擦单元和黏弹性单元并联组成，如图4所示。

图4中Fe，Ff，Fv和F分别为弹性单元、摩擦单元、黏弹性单元和衬套总模型的响应力，N。

2.1 弹性单元和摩擦单元

衬套轴向弹性特性在线变形范围内符合胡克定律［10］：

式中：Fe为弹性力，N；x为衬套变形量，mm；Ke为静态弹性刚度，N·mm-1。在振幅为 x0的正弦激励下，弹性力的幅值Fe0为Kex0，且无能量损耗。

图4 橡胶衬套动力学模型组成示意图

摩擦单元中摩擦力与位移的关系为［10］

式中：Ff为摩擦力，N；x为加载位移，mm；Ffmax为最大摩擦力，N；x2为Ff从0开始，增加至 Ffmax／2时的位移，mm；（xs，Ffs）为摩擦力随位移变化曲线上的状态参考点，此处取（0，0）；sign（x·）为位移的变化方向，位移增加时为正，反之为负。在振幅为x0的正弦激励下，摩擦力的幅值和损耗的能量分别为

式中：Ff0为摩擦力的幅值，N；Ef为一个周期损耗的能量，N·mm；a0＝Ff／Ffmax。

2.2 黏弹性单元

基于传统黏弹力模型，Gemant等［11］学者提出了分数导数模型，因能反映加载历程对黏弹性材料动态特性的影响，得到广泛应用。最常见的分数导数Kelvin-Voigt模型和Maxwell模型，进一步提出的分数导数Zener模型以及在此基础上提出的包含两个不同分数阶的分数导数模型，均为最高阶次不大于1的低阶分数导数模型，在模拟橡胶材料随频率的变化趋势中的损耗模量存在一定误差［6］。高阶（模型参数取不同值时最高阶次大于1）分数导数模型对黏弹性力学特性的描述更加精确［8］，因此本文中通过将分数导数Kelvin-Voigt模型和Maxwell模型并联，建立高阶分数导数模型描述橡胶材料的黏弹特性，如图5所示。

图5 高阶分数导数黏弹力模型示意图

图5 中，k1和 c1分别为分数阶Kelvin-Voigt模型中弹簧的弹性模量和黏壶的黏性系数，k2和c2分别为分数阶Maxwell模型中弹簧的弹性模量和黏壶的黏性系数。k1和k2反映了材料的弹性，c1和c2反映了材料的黏性。

根据图5，黏弹力与位移的关系为

式中：λ1＝k1／c1；λ2＝k2／c2；Fv（t）为黏弹性力，N；x（t）为加载位移，mm；α，β和 γ为分数导数阶数，取值范围均为（0，1），其中最高阶数为 α+γ（最大值为2）。对式（5）进行傅里叶变换：

则该分数导数黏弹力单元确定的复刚度为

式中：K*v（ω）为黏弹力单元的复刚度，N·mm-1；将

进一步推导出，在振幅为x0的正弦激励中，响应力实部和虚部的幅值Fv0Re和Fv0Im分别为

一个加载循环中的能量损耗为

式中Ev为一个加载周期的能量损耗，N·mm。

2.3 橡胶衬套动力学模型

将弹性单元、摩擦单元与黏弹性单元并联，得到衬套动力学模型响应力和加载位移的关系为

式中：F（x）为衬套模型的响应力，N；Fe（x），Ff（x）和Fv（x）分别为弹性单元、摩擦单元和黏弹性单元的响应力，N；x为加载位移，mm。

衬套模型响应力振幅为式中：Fe0和Ff0分别为弹性单元和摩擦单元响应力振幅，N；Fv0Re和Fv0Im分别为黏弹性单元响应力实部和虚部振幅，N。

一个加载循环中，衬套模型的能量损耗为

式中：Ef和Ev分别为一个加载循环中摩擦单元和黏弹性单元的能量损耗，N·mm。

振幅为x0的正弦加载工况下，衬套的动刚度和阻尼系数可通过式（15）和式（16）计算得到：

式中：Kdyn为衬套模型的动刚度，N·mm-1；D为衬套模型的阻尼系数，（N·s）·mm-1。

3 橡胶衬套力学模型参数识别

3.1 弹性单元和摩擦单元参数的识别

选取频率为1 Hz、振幅为2 mm正弦加载试验工况作为两个单元参数的识别工况，该工况下响应力F与位移x的迟滞环曲线如图6所示。

图6 低频加载工况下位移与力的关系曲线

根据图6，位移接近极限位置处曲线的斜率可近似表示该衬套弹性单元的静态弹性刚度Ke。延长接近极限位置的上、下两条切线，摩擦力模型中最大摩擦力Ffmax的值即为两条切线间的垂直距离的1／2。取曲线的最大斜率的值Kmax，由式（17）可求得摩擦力单元中的x［10］。

弹性单元和摩擦单元参数识别结果见表1。

表1 弹性单元和摩擦单元参数识别结果

3.2 黏弹性单元参数的识别

黏弹性单元参数 k1，c1，k2，c2，α，β和 γ7个参数通过小振幅激励下衬套动态特性获得［10］，为充分体现衬套在不同频率下的非线性黏弹特性，选取振幅为 0.1 mm、频率分别为 1，2，4，8，12，16，20，30和40 Hz的9个正弦加载试验工况进行7个参数的识别。

黏弹性单元参数识别实质上为最优参数估计问题，为使识别结果尽可能接近试验数据，设置衬套动力学模型计算的衬套动态参数值与试验值的误差最小的优化函数来实现。因此，建立基于动刚度和阻尼系数误差优化的目标函数：

为保证黏弹性单元参数识别的准确性，提高橡胶衬套模型的整体预测精度，参数识别过程中，需确保各频率工况下模型计算得到的与试验值均不能有较大误差。一般认为模型的计算精度在90%以上，即可满足工程应用需要，因此，建立约束条件：

由于高阶分数导数黏弹性单元为复杂的非线性模型，且需要识别的参数个数较多，粒子群优化（particle swarm optimization，PSO）算法，作为一种非解析寻优算法，具有大范围收敛性（不需要较好的近似值作为搜索的起点，有效避免对初值的敏感性）、高效并行、鲁棒性和通用性强等优点，对非线性系统参数识别效果较好，且能够实现多参数识别［12-13］，因此，把PSO算法引入到高阶分数导数黏弹力模型参数识别中。为提高算法性能，避免陷入局部极值点，提出一种基于动态邻居广义学习策略的粒子群优化（dynamic neighbor and generalized learning particle swarm optimization，DNGL-PSO）算法的高阶分数导数黏弹性单元参数识别方法。通过在PSO算法中引入邻居动态组建策略和广义学习策略，来降低参数识别过程中算法陷入局部最优解的可能性，使算法有较大的概率向全局最优解收敛，从而提高黏弹性单元参数的识别精度。

在DNGL-PSO算法运行过程中每一迭代时刻，各粒子根据自身适应度函数值动态地组建邻居，选择与其欧氏距离较近的粒子作为邻居，邻居个数为群规模的 1／4～1／3［14］。同时，引入一种广义学习策略，使各粒子“个体学习”部分的学习样本为其所有邻居（包含其自身历史最优位置），保证粒子的每一维学习对象在当前迭代时刻为最优。粒子i第d维的学习样本pdbn_i的选取规则为

式中：pi和Fit（pi）分别为粒子 i的历史最优位置及其所对应的适应度值；n为种群规模；Ri为粒子i的邻居构成的集合。

粒子i的速度和位置更新公式为

式中：t为当前迭代次数；v-i（t）为粒子i在第t代的速度；x-i（t）粒子i在第t代的位置；ω为粒子的惯性权重；ε1和ε2为粒子的学习因子；r1和r2为相互独立在［0，1］均匀分布的随机数；pg（t）为整个粒子群迄今为止搜索到的最优位置；pbn_i（t）为“个体学习”部分粒子i每一维的学习样本。其中，ω为线性递减。

以式（18）所示目标函数为适应度函数Fit（t），采用DNGL-PSO算法对黏弹性单元参数识别，经过多次调试，最终算法参数设置如下：取种群规模n为48，每个粒子的邻居个数为16；惯性权重ω为0.9～0.4线性递减；学习因子 ε1和 ε2均为2；r1和 r2为在［0，1］均匀分布的随机数；粒子维数d为7，位置向量xi＝（k1，c1，k2，c2，α，β，γ），考虑所建高阶分数导数单元的特点和参数识别结果的准确性，设置粒子各维搜索空间下限xmin＝（0，0，0，0，0，0，0），上限xmax＝（103，103，103，103，1，1，1）；最大迭代次数 Gmax取200。

基于DNGL-PSO的高阶分数导数黏弹性单元参数识别流程如图7所示，识别结果如表2所示。

表2 黏弹性单元参数识别结果

4 橡胶衬套动力学模型验证及修正

4.1 橡胶衬套动力学模型验证

模型参数识别完成后，基于衬套动力学模型进行了不同频率和振幅的正弦激励动态仿真计算，结果如图8和图9所示，当前模型预测的衬套动刚度和阻尼系数与频率和振幅的关系均与试验结果（图2和图3）一致。由于摩擦单元的能量损耗不受频率影响，既使在激励频率接近于零时，由模型计算出的衬套能量损耗仍呈现一定的值。因此，该模型能够较为准确地反映橡胶衬套轴向动态特性分布特征。

图7 基于DNGL-PSO的黏弹性单元参数识别流程图

图8 衬套动刚度预测结果

预测值与试验结果的对比如图10和图11所示，基于小振幅试验工况识别出的黏弹性单元参数，使所建立的高阶分数导数黏弹性衬套模型能够在较宽频率范围内精确地预测衬套的轴向动态特性。但对于大振幅加载试验工况，预测精度有所下降。

4.2 橡胶衬套动力学模型修正与验证

试验和研究表明，除频率外，振幅对衬套的黏弹特性也有一定影响［2］。为提高模型对大振幅试验工况的预测精度，对当前模型进行修正。分数导数黏弹性模型中，分数导数阶数是描述黏性强弱的物理参数，对于特定材料，阶数是固定的［15］。因此，进一步地，将黏弹性单元的分数导数系数k1，c1，k2和c2设为与加载振幅相关，对当前模型进行修正。

令小振幅工况下识别出的分数导数阶数α，β和γ保持不变，采用DNGL-PSO算法，根据振幅为0.1，0.8和2 mm的不同加载频率（每个振幅下均包含1，2，4，8，12，16，20，30和 40 Hz 9个频率工况）的试验结果对单元中的k1，c1，k2和c2重新进行识别，并使用二次多项式拟合各参数与振幅x0的关系：

基于新修正的高阶分数导数黏弹性并联衬套动力学模型，进行正弦激励动态仿真计算，与试验结果进行对比，如图12和图13所示。可见最终建立的黏弹性高阶分数导数并联动力学模型在各种试验工况下均具有较高的预测精度。

图12 衬套动刚度模型预测值与试验值对比

图13 衬套阻尼系数模型预测值与试验值对比

在振幅为0.1～2 mm之间，频率在1～40 Hz之间，分别对初期模型、最终模型预测结果与试验结果进行统计，得到两个模型对各加载振幅工况的衬套动刚度和阻尼系数最大预测误差的平均值和最大值，如表3所示。

由表3可看出，最终模型能够准确描述预测工况内橡胶衬套轴向的动态特性，其中，动刚度的平均误差和最大误差分别为1.99%和2.82%，阻尼系数的平均误差和最大误差分别为5.03%和5.76%。

表3 初期模型和最终模型对各加载振幅下最大预测误差的对比结果 %

文献［3］中包含低阶分数导数黏弹性单元的衬套模型，对阻尼系数预测最大误差为10.43%，本文中所建模型显著降低了模型在预测工况内的波动，说明高阶分数导数黏弹性单元能够有效减小模拟橡胶材料随频率变化的阻尼特性的误差。文献［16］中基于多个黏弹性单元叠加的衬套模型，对动刚度和阻尼系数最大预测误差分别为6%和8%，精度较高，但对衬套材料物理本质的描述不够准确，本文中所建衬套模型具有很高的整体预测精度和适用性，能更加接近地描述衬套材料的物理本质。