一种磁流变弹性体模型参数识别新方法及其应用研究

王 鹏， 杨绍普， 刘永强， 董旭峰， 赵义伟

(1.石家庄铁道大学 省部共建交通工程结构力学行为与系统安全国家重点实验室，石家庄 050043；2.石家庄铁道大学 交通运输学院，石家庄 050043； 3.大连理工大学 材料科学与工程学院，大连 116024)

磁流变弹性体(magnetorheological elastomers，MRE)是一种磁敏智能材料，由弹性基体、磁性颗粒和添加剂组成[1]，具有良好的磁流变效应。与磁流变液相比，磁流变弹性体由于没有颗粒沉降而具有稳定的力学和磁流变性能，在调谐阻尼器、节点衬套、变刚度悬架等领域具有广泛的应用前景[2-3]。

目前，对磁流变弹性体力学性能的研究主要从微观和宏观两个方面进行。基于磁偶极子理论提出了微本构模型，主要描述了磁流变效应的机理以及磁致伸缩模量与各种因素之间的关系。Jolly等[4]建立了磁偶极子模型，该模型假设磁相互作用仅存在于链内相邻粒子之间。Davis[5]基于有限元分析计算剪切模量，假设颗粒均匀地排列成一条长链，提出了链状模型。但链状模型忽略了颗粒分布不均匀性对磁流变效应的影响。此后，一些学者又提出了三维结构模型，如柱模型和网格模型。Chen等[6]基于MRE微观形貌扫描照片的统计数据建立了有限柱模型。Zhang等[7]针对两种特定的颗粒分布模式(点阵型和体心型)提出了一种网格计算模型。三维结构模型虽然能够精细地反映磁性颗粒在MRE中的分布结构，但通常构造复杂。

为使磁流变弹性体应用于工程实际，还需要可以描述其宏观力学行为的力学模型。宏观力学模型应综合考虑外加磁场、激励幅值以及加载频率等对MRE力学性能的影响。Li等[8]提出了四参数线性黏弹性模型预测MRE的力学特性，但该模型不能有效地描述MRE的非线性滞后特性。Kari等[9]认为大应变下MRE滞回曲线的平行四边形特征是由磁性颗粒与基体之间的摩擦引起的，故采用库伦摩擦模型描述此现象。然而，包含库仑摩擦单元的模型不能准确反映加载过程中切线模量的逐渐变化规律。作为描述磁滞回线的经典模型，Bouc-Wen模型已被广泛用于描述与磁场相关的非线性磁滞系统。Dominguez等[10]将磁流变阻尼器的励磁电流引入Bouc-Wen模型，并考虑了电流变化对磁流变阻尼器机械性能的影响。Yang等[11]建立了Bouc-Wen模型来描述MRE隔离器的机械特性，该模型能够准确描述MRE隔离器的滞后特性。赵义伟等[12]针对非识别激励工况下模型预测的阻尼力误差较大问题，提出了一种描述滞回特性的改进Bouc-Wen模型。Wang等[13]考虑了MRE的非线性、频率依赖性，提出了基于Bouc-Wen的磁流变弹性体非线性本构模型。上述Bouc-Wen模型能够准确模拟MRE的力学行为。

Bouc-Wen模型的参数辨识方法很多，可分为两大类。一类是通过分析Bouc-Wen模型的强迫极限环特征来识别参数。Wang使用Bouc-Wen模型再现了隔振器的迟滞力学特性，并利用极限环方法识别了模型参数[14]。Rodriguez等[15]用Bouc-Wen模型研究磁流变阻尼器动态力学性能，通过分析强迫极限环特征对模型参数进行识别。极限环方法使用多个单独的数据点来识别模型参数，这使得结果受到数据噪声的影响，并且对试验采样频率有很高的要求。另一种是采用具有适当适应度函数的数值优化算法。该方法采用不同的数值优化算法来识别多个参数。刘永强等[16]采用遗传算法对磁流变阻尼器Bouc-Wen模型进行参数辨识，并采用缩小参数取值范围的方法逐渐提高遗传算法的识别精度。Giuclea等[17]提出了一种改进的Bouc-Wen动力学模型，通过遗传算法获得模型参数。随后，Charalampakis等[18]采用混合进化算法对Bouc-Wen滞回系统进行辨识。由于数值优化算法的应用，结果很容易陷入局部最优解。此外，该方法忽略了模型参数的物理意义。

在本研究中，采用Bouc-Wen模型定量地描述MRE的力学特性。分析了模型参数与应变幅值、频率和磁场之间的耦合关系，提出了多工况参数识别方法，得到了Bouc-Wen模型的参数表达式。最后，利用识别和非识别工况下MRE的试验数据对模型识别结果进行验证，证明了多工况参数识别方法的有效性。

1 MRE制备及性能测试

选用天然橡胶(natural rubber，NR)、羰基铁粉(carbonyl iron particles，CIP)和添加剂作为原材料，通过塑炼、混合和硫化等工艺制备了CIP质量分数80%的各向同性MRE样品。CIP平均粒径5.5 μm，使用的添加剂主要有硬脂酸(Stearic Acid，SA)、氧化锌(ZnO)、促进剂(CZ)、防老剂(RD)、防老剂(4010NA)、古马隆树脂和硫磺(S)，MRE成分配比，如表1所示。

表1 MRE成分配比Tab.1 The composition ratio of MRE sample

采用Physica MCR301流变仪(Anton Paar，奥地利)对MRE进行动态性能试验。测试时，首先将试件放置于下平台的中心，随后上部转子向下运动直至与试件上表面接触，并施加法向力(3 N)使试件与上下极板紧密贴合，防止试件发生滑动。旋转流变仪通过改变线圈电流实现在平行板系统上施加不同的磁场，上部转子按照设定程序施加平行于试件的旋转剪切应变激励，并由内部精密传感器测定试件反馈的应力，经过仪器控制软件计算得到所需结果。在振荡剪切测试模式下，每个正弦应变加载周期内软件可记录257个应力应变数据点，根据这些数据点绘制每个正弦应变加载周期内的应力应变曲线。在正弦剪切应变激励下测试了MRE在不同应变幅值、频率和磁场条件下的稳态响应，试验工况如表2所示。图1为旋转流变仪、试件及测试示意图，其中试件为直径20 mm，厚度1 mm的MRE圆片。

表2 MRE测试工况Tab.2 MRE test conditions

图1 旋转流变仪、试件及测试示意图Fig.1 Rotary rheometer, specimen and schematic of the measuring system

2 MRE力学模型及其验证

2.1 Bouc-Wen模型

MRE具有明显的非线性滞后特性，线性黏弹性参数模型不能模拟MRE磁滞回线的平行四边形特性。作为描述磁滞回线的经典模型，Bouc-Wen模型能够很好地模拟这一特性，特别是在大应变幅值下。Yang等首先利用Bouc-Wen模型拟合MRE隔离器的力-位移关系，该模型由Bouc-Wen滞回算子、弹簧和阻尼元件并联组成，如图2所示。Bouc-Wen模型结构简单，通用性强，且能够很好地模拟滞回特性。本文采用Bouc-Wen模型来模拟MRE的力学特性。控制方程为

(1)

(2)

图2 Bouc-Wen模型Fig.2 The modified Bouc-Wen model

式中：k表示弹性系数;c表示黏性系数;参数z为中间变量;α介于0～1之间，表征模型的线性程度;A、n、β、γ是无量纲参数，共同决定了滞回曲线的形状和大小。

由于Bouc-Wen模型存在微分项，采用传统方法数值仿真较困难，根据控制方程在Simulink软件中建立了仿真模型。

2.2 GA-PSO混合算法

粒子群优化算法(particle swarm optimization，PSO)是假设在D维的搜索空间中，有M个粒子构成一个群体，当前的迭代次数记为λ，粒子当前位置记为Xi，当前飞行速度记为Vi，粒子当前最优位置记为qi，粒子群当前最优位置记为g。每个粒子的位置就是一个潜在的解，代入目标函数L计算适应值，根据适应值大小衡量解的优劣。第i个粒子最优位置更新

qi(λ+1)=

(3)

粒子群当前最优位置由下式确定

L[g(λ)]=min{L[q1(λ)],…,L[qm(λ)]}

(4)

粒子在搜索空间中的速度更新公式

vij(λ+1)=ηvij(λ)+d1r1[qij(λ)-xij(λ)]+

d2r2[gj(λ)-xij(λ)]

(5)

粒子在搜索空间中的位置更新公式

xij(λ+1)=xij(λ)+vij(λ+1)

(6)

其中，i=[1,2,…,M]，j=[1,2,…,D]；η表示惯性权重；d1和d2是加速常数；r1和r2为0～1内的随机数。设定粒子的速度区间(vmin，vmax)，位置范围(xmin，xmax)，如果超出范围则取边界值。

PSO算法具有较好的优化能力，但局部搜索性能较差，导致后期收敛速度减慢，产生不准确的结果。遗传算法(genetic algorithm，GA)作为一种基于概率的优化方法，能够自适应地调整搜索方向，具有出色的局部搜索能力，但优化过程耗时且对种群初始化敏感。GA和PSO算法结合使用，既改善了GA算法易陷入局部极值、运算速度慢的缺点，又提高了PSO算法的运算精度[19]。因此，采用GA-PSO混合算法对Bouc-Wen模型参数进行识别。首先初始化粒子群种群，然后计算每个粒子的适应度值。根据粒子适应度值更新个体最优位置和群体最优位置。在PSO迭代产生新一代个体后，再通过遗传操作(交叉和变异)将当前种群进行更新[20]。设置惯性权重η=1，加速常数d1=d2=3.5，粒子数M=36，迭代次数400，染色体编码长度7，交叉概率0.7，变异概率0.3。利用Liu等[21]提出的误差计算公式作为参数识别的目标函数L，公式如下：

(7)

MRE力学模型的参数识别是一个最小化目标函数的优化过程。流程图如图3所示。

图3 基于GA-PSO算法的参数识别流程Fig.3 Flowchart of Hybrid GA-PSO algorithm

2.3 MRE力学模型验证

分析Bouc-Wen模型能否准确表征不同工况下MRE的力学特性。利用GA-PSO优化算法和MRE试验数据，按照识别流程进行参数辨识。

(1) 磁感应强度为0，应变频率为1 Hz，不同应变幅值下，MRE的仿真与试验结果对比，如图4所示，吻合率如表3所示。

(a)

(b)图4 不同应变幅值下的仿真与试验结果对比Fig.4 Simulation data and test data under different amplitudes

由图4可知，随着应变幅值增加，MRE的滞回环形状由倾斜的椭圆逐渐演变成两头尖的橄榄形，进一步发展呈现出平行四边形特征，而整个变化过程中仿真曲线与试验曲线均吻合较好。需要指出的是，当应变幅值为15.80%和29.3%时，试验曲线有部分重叠，而仿真曲线未出现重叠，这一现象是因为MRE的储能模量会随着应变幅值增加会快速下降，所以试验曲线出现了重叠；对于仿真曲线而言，拟合精度与选取的MRE力学模型有关，也和参数识别结果有关，所以仿真得到的滞回曲线并未出现重叠现象，但是整体变化趋势与试验曲线保持一致。说明Bouc-Wen模型可以准确描述MRE非线性滞回特性，反映黏弹性随应变幅值的变化规律。

(2) 磁感应强度为405 mT，应变幅值为0.40%，不同应变频率下，MRE的仿真与试验结果对比，如图5(a)所示，吻合率见表3。

由图5(a)可知，磁流变弹性体的动态黏弹性对加载频率有一定的依赖性。随着应变频率增加，试验数据绘制的滞回环的斜率和面积均不断变大。四种应变频率下仿真数据很好地反映了这一力学特性，且与试验结果吻合良好。说明Bouc-Wen模型能有效追踪滞回环随频率的变化规律。

(a) 不同频率

(b) 不同磁感应强度图5 仿真结果与试验结果对比Fig.5 Comparison of the numerical and test results

(3) 应变幅值为2.51%，应变频率为1 Hz时，不同磁感应强度下，MRE的仿真与试验结果对比，如图5(b)所示，吻合率见表3。

表3 Bouc-Wen模型参数识别结果及吻合率Tab.3 Identification values of the Bouc-Wen model parameters and fitness values

由图5(b)可知，随着磁感应强度的增加，试验数据的滞回环包围的面积与对角线斜率逐渐变大，即MRE的动态黏弹性具有磁场依赖性，不同磁场下仿真结果与试验结果吻合效果较好，二者滞回曲线差异不大。此外，对于试验值，磁感应强度由456 mT增大到545 mT时滞回环的最大应力值的增长幅度明显小于磁场强度增长初期时滞回环最大应力值的增长幅度，这是MRE磁致模量增加接近磁场饱和引起的；图中可以看出仿真数据的滞回曲线同样准确地反映了MRE磁致黏弹性的磁场饱和特性。说明Bouc-Wen模型可以准确反映磁流变效应对MRE力学性能的影响。

3 多工况参数识别方法

Bouc-Wen模型共有7个参数，若每次应用模型时均要对7个目标参数寻找最优结果，不仅影响模型使用的便捷性，且很难保证优化结果的一致性。因此，需要提供确定MRE力学模型参数的识别方法。

3.1 参数耦合关系分析

(8)

(9)

将式(8)代入式(9)可得：

(10)

进一步分析，当应变频率和磁场相同时，将应变幅值1.36%工况得到的参数识别结果用于预测其他应变幅值(2.51%、4.64%、8.57%、15.8%、29.3%)工况，发现仿真值与试验值拟合误差较大，即单一应变幅值工况得出的识别结果不能准确预测其他幅值工况的滞回环。同样，将其他应变幅值(2.51%、4.64%、8.57%、15.8%、29.3%)工况的参数识别结果用于预测应变幅值1.36%工况时，仿真值与试验值吻合度也较低，说明Bouc-Wen模型参数与应变幅值存在耦合。

图6(a)为磁感应强度0，应变频率1 Hz时，不同应变幅值下的MRE试验数据的滞回曲线。由图可知，随着应变幅值增加，滞回曲线的等效刚度(滞回曲线对角线斜率)不断减小，滞回曲线的等效阻尼(滞回曲线的面积来表征[22])不断增大，而参数k表示弹性系数，c表示黏性系数，也从另一侧面反映出参数k、c与应变幅值相关。

图6(b)表示应变频率为1 Hz，应变幅值为2.51%时，不同磁感应强度下MRE试验数据滞回曲线。由图可知，不同磁感应强度下，滞回曲线的形状基本保持不变，可近似视为滞回曲线形状不随磁场改变而改变，而滞回曲线的形状由参数A、γ、β、n共同调控，所以参数A、γ、β、n受磁场影响较小，此处不再考虑其与磁场之间的耦合关系。此外，随着磁感应强度的增加，滞回曲线的等效刚度和等效阻尼不断增大，说明参数k、c与磁场存在耦合。参数α影响滞回曲线的线性程度，与滞回曲线的刚度、阻尼相关，因此与磁场存在一定的函数关系。

(a) 不同应变幅值

(b) 不同磁感应强度图6 MRE应力应变的试验数据Fig.6 Test results of MRE stress and strain

综上分析，Bouc-Wen模型参数与频率、幅值及磁场存在耦合。与频率相关的参数有c；与幅值相关的参数有k、c、α、A、γ、β、n；与磁场相关的参数有k、c、α。

3.2 模型参数识别过程

上节分析可知，Bouc-Wen模型参数同时与多个因素存在耦合，要将所有的耦合关系量化表示难度较大。为了降低参数识别难度，本文对参数耦合关系进行适当简化，提出了多工况参数识别方法。

基本思想是：在整个识别过程中忽略磁场和应变频率对α的影响；同时先不考虑应变频率对k、c的影响，即在任一应变频率下探寻模型参数与幅值、磁场的函数关系；再结合应变频率对k、c的影响，对其参数表达式进行适当修正，以提高不同应变频率下的拟合精度。具体步骤如下：

(1) 应变频率为1 Hz时，对25种工况(应变幅值取1.36%、2.51%、4.64%、8.57%、15.8%；磁场取0、115 mT、234 mT、456 mT、545 mT)的试验数据进行第一次参数识别，得到kij、cij、αij、Aij、γij、βij、nij；i表示幅值，j表示磁场。

(6) 结合频率f对参数k和c的影响，对参数k和c的表达式进行修正，得到[k,c]=F3(f)·F2(xm,B)。

根据上述过程，利用步骤(1)中25种工况的试验数据对Bouc-Wen模型进行参数识别，分别得到第一次识别结果见表4，第二次识别结果见表5。此处需要指出的是在B=0，f=1 Hz时，应变幅值变化下表3和表4中数据是先后进行两次不同的识别得到的结果，可以看出数据的吻合率相差不大，但数值不同，这也从侧面印证了不同的识别过程尽管拟合误差相近，但很难保证优化结果一致。因此，研究MRE力学模型参数的识别方法是非常必要的。

表4 Bouc-Wen模型参数第一次识别结果Tab.4 First identification values of the Bouc-Wen model parameters

表5 Bouc-Wen模型参数第二次识别结果Tab.5 Second identification values of the Bouc-Wen model parameters

根据表4中不同应变幅值下参数α、A、γ、β、n的均值，利用Matlab软件的曲线拟合工具箱确定参数α、A、γ、β、n与应变幅值xm的函数关系。参数表达式如式(11)所示，拟合结果如图7所示。

0.184 8xm+0.713 9

12.23xm+38.73

0.243 2xm+0.826 9

(11)

(a)

(b)

(c)

(d)

(e)图7 Bouc-Wen模型参数拟合结果Fig.7 Fitting results of Bouc-Wen model parameters

根据表5，确定参数k、c与应变幅值xm和磁感应强度B的函数关系。参数表达式如式(12)所示，拟合结果如图8所示。

9.898×10-9xmB3-1.718×10-11B4

(12)

(a)

(b)图8 Bouc-Wen模型参数k、c拟合结果Fig.8 Fitting results of Bouc-Wen model parameters k, c

选取应变幅值2.51%，利用式(11)计算参数α、A、γ、β、n，得到的结果带入Bouc-Wen模型，使模型只剩余未知参数k、c。再利用磁感应强度234 mT，应变幅值2.51%，频率分别为0.5 Hz、1 Hz、3 Hz、5 Hz的试验数据对参数k、c进行二次识别，表6为不同频率下参数k、c的识别结果(各频率工况的绝对值除以1 Hz工况的绝对值即为相对值)。

表6 不同频率下参数k和c识别结果Tab.6 Identification values of the parameters k and c under different frequencies

由表6中参数k、c的相对值，确定参数k、c关于应变频率f的函数表达式

Fk(f)=-0.246 5f3+2.077f2-4.465f+3.634

Fc(f)=-0.336 9f3+3.046f2-8f+6.29

(13)

再利用式(13)对式(12)进行修正，得参数k、c关于应变幅值xm、应变频率f、磁感应强度B的函数表达式为

k(xm,B,f)=Fk(f)×(9.687-4.779xm-3.034×

5.004×10-6xmB2-1.064×10-7B3+2.748×

5.354×10-9xmB3)

c(xm,B,f)=Fc(f)×(7.073-1.305xm-5.492×

9.197×10-6xmB2-1.249×10-7B3+1.854×

1.718×10-11B4)

(14)

4 参数识别结果验证

分析Bouc-Wen模型的仿真结果与试验结果的吻合情况，验证多工况参数识别方法所得结果的正确性。

4.1 识别工况下验证

对应变频率1 Hz，不同应变幅值(1.36%、2.51%、4.64%、8.57%、15.8%)和不同磁场(0、115 mT、234 mT、456 mT、545 mT)的25种识别工况，仿真与试验结果的吻合率见表7，曲线对比见图9。

表7 识别工况下仿真和试验结果的吻合率Tab.7 Fitness values of the numerical and test results

(a)

(b)

(c)

(d)

(e)图9 识别工况下仿真和试验结果对比(f=1 Hz)Fig.9 Comparison of the numerical and test results under identification conditions (f=1 Hz)

由图11和表7可知，无论是不同幅值，还是不同磁感应强度下，仿真结果与试验结果的吻合率都高于93.2%，模型仿真的滞回曲线能准确模拟MRE的力学特性，说明参数识别结果对识别工况试验数据有效。

4.2 非识别工况下验证

利用非识别工况的MRE试验数据对Bouc-Wen模型参数识别结果进一步验证。

对于磁感应强度234 mT时，应变幅值1.36%、4.64%，应变频率分别取0.5 Hz、3 Hz、5 Hz的试验工况；以及磁感应强度115 mT时，应变幅值取1.36%、2.51%、4.64%，应变频率分别为0.5 Hz、3 Hz、5 Hz的试验工况，仿真与试验结果的吻合率见表8，曲线对比见图10。

表8 非识别工况下仿真和试验结果的吻合率Tab.8 Fitness values of the numerical and test results

由图10和表8可知，不同应变幅值、频率和磁感应强度作用下，仿真结果与试验结果的滞回曲线吻合良好，说明预测模型能够准确的模拟MRE滞回特性。由表8数据可知，各工况的吻合率均保持在94%以上，说明参数识别结果正确，证明了本文提出的多工况参数识别方法的有效性。

(a) 磁场B=234 mT，幅值xm=1.36%

(b) 磁场B=234 mT，幅值xm=4.64%

(c) 磁场B=115 mT，幅值xm=2.51%

(d) 磁场B=115 mT，幅值xm=1.36%

(e) 磁场B=115 mT，幅值xm=4.64%图10 非识别工况下仿真和试验结果对比Fig.10 Comparison of simulation and test results under non identification conditions

5 结 论

本文采用工业炼胶方法制备了CIP质量分数80%各向同性的天然橡胶基MRE，利用流变仪在正弦剪切应变激励对MRE进行了测试。采用Bouc-Wen模型来表征MRE的非线性滞回特性，试验数据验证了模型的有效性。分析了Bouc-Wen模型参数与应变幅值、频率和磁场的耦合关系分析，提出了多工况参数识别方法，将仿真值与试验值之间的误差作为目标函数，利用GA-PSO混合算法对模型参数进行识别，采用多种应变幅值、频率和磁场下的试验数据对识别结果进行验证。结论如下：

(1) Bouc-Wen模型能够准有效拟不同应变幅值、频率和磁场下磁流变弹性体的力学特性。在磁感应强度0～545 mT，应变幅值0.04%～15.8%，应变频率 0.5～5 Hz的加载工况，模型的拟合精度较高，可以准确描述MRE非线性滞回特性，反映MRE的动态黏弹性随着应变幅值的变化规律，追踪滞回环随着频率的变化趋势，以及磁致黏弹性的磁场饱和特性。

(2) 提出了多工况参数识别方法，得到了Bouc-Wen模型参数表达式，量化了参数与影响因素之间的关系。识别和非识别工况下Bouc-Wen模型仿真结果与试验结果的吻合率均高于93%，证明了多工况参数识别方法的有效性，该方法不限于所使用的Bouc-Wen模型，可广泛用于MRE其他力学模型的参数识别。