磁流变阻尼器具有体积小、能耗低、结构简单、阻尼力大、动态范围广、频响高、适应面大等特点，且能够依据系统振动实时调节阻尼力、改善系统振动响应，因而在振动控制领域具有广阔的应用前景.

准确的磁流变阻尼器出力预估模型对于磁流变阻尼器的优化设计具有重要意义. 由于磁流变效应的复杂性、磁流变液本身的非线性及非均匀磁场的影响，磁流变阻尼器的出力难以准确地估计. 现在通常用Bingham黏塑性模型^[1]描述磁流变液的应力-应变关系，但由于磁流变液在屈服前表现为黏弹性，国内外的学者引入黏弹性流体本构对磁流变液进行研究. Ferrás等^[2-3]基于Maxwell黏弹性流体本构模型，通过分数阶微积分建立磁流变液的黏弹性本构关系. Peng等^[4]基于微观流体动力学，建立磁流变液的黏弹性模型. 蔡路^[5]基于Oldroyd-B黏弹性本构仿真，分析阻尼器阻尼通道内磁流变液流动特性. 为了模拟非均匀磁场下磁流变液的黏度分布，许多学者采用有限元多场耦合仿真的方法，将磁场仿真得到的非均匀磁场导入流场中，根据磁流变液的材料属性得到此时的流体表观黏度分布，通过流固耦合分析得到活塞运动时的阻尼器出力. Paul等^[6]通过多场耦合仿真分析双黏性本构下的磁流变阻尼器出力. Kazakov等^[7]根据双黏性本构对磁流变阻尼器的磁场、流场与热场进行耦合分析，得到阻尼器活塞运动时受到的阻力. Sadak等^[8]基于Bingham黏塑性本构模型，通过磁场仿真与理论计算分析不同电流及活塞速度下的阻尼器出力情况. Xin等^[9]通过磁热仿真与理论公式，分析不同电流与活塞激励下的阻尼器出力. 于振环等^[10]通过磁-流和流-固耦合仿真，分析阻尼器在非控及通电情况下的流场压力及磁场、示功特性.

基于阻尼器出力模型，可以分析阻尼器各结构参数对阻尼器出力的影响，对磁流变阻尼器进行优化设计. Fathima等^[11]基于利用磁流变液的Bingham本构模型，利用有限元软件分析阻尼间隙对阻尼器出力及可调范围的影响. Gurubasavaraju等^[12]通过耦合场有限元仿真分析Bingham流体本构下，阻尼间隙、激励频率与活塞速度对阻尼器出力的影响. Parlak等^[13]以磁流变阻尼器的目标阻尼力和最大磁通密度为目标，通过有限元多场耦合进行设计优化.

现在，磁流变阻尼器的有限元出力仿真大多是基于简单的黏塑本构模型建立的. 对于黏弹性模型，现在大都应用于阻尼器动力学参数建模阶段，这是基于阻尼器的示功实验数据进行辨识参数建立的数学模型，不能用于阻尼器的设计初期. 为了在阻尼器设计初期建立准确的阻尼器出力模型，本文参考阻尼器动力学模型^[14-15]，以Bingham黏塑性模型与Maxwell（Kelvin）黏弹性模型为基础，在Bingham模型的基础上串联（并联）一个线性弹簧单元建立串联（并联）黏弹塑性本构模型. 根据建立的磁流变液本构模型，通过COMSOL软件建立阻尼器的出力模型. 基于该模型分析各结构参数对磁流变阻尼器力学特性的影响，计算各参数对阻尼器出力的灵敏度，为磁流变阻尼器的结构优化提供了理论依据.

Bingham黏塑性模型由理想刚塑性模型和牛顿黏壶并联而成，如图1所示，本构关系^[16]为

(1) $ \left. {\begin{array}{*{20}{l}} {\tau = \tau (H){\rm{sgn}}{\kern 1pt}{\kern 1pt}(\dot \gamma ) + \eta \dot \gamma ,{\kern 1pt} {\kern 1pt} {\kern 1pt}{\kern 1pt}{\kern 1pt} \;|\tau |{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \text{＞} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} |{\tau _0}|;\;}\\ {\dot \gamma = 0,{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} |\tau |{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \text{＜} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} |{\tau _0}|{\rm{.}}} \end{array}} \right\} $

式中：τ₀为磁致剪切屈服强度， $\dot \gamma $为流体的切应变率，H为磁场强度，η为与磁场无关的液体屈服后黏度. 应力-应变关系如图1所示.

为了描述磁流变液的黏弹塑性，在Bingham黏塑性模型的基础上增加一个线性弹簧单元. 下面介绍2种基于Bingham模型的黏弹塑性模型.

1）并联黏弹塑性模型（并联模型）. 并联黏弹塑性模型由一个线性弹簧与一个Bingham黏塑性模型并联而成，如图2所示，本构关系为

(2) $\left. {\begin{array}{*{20}{l}} {\gamma {\kern 1pt} {\kern 1pt}{\rm{ = }}{\kern 1pt} {\kern 1pt}{\gamma _1}{\kern 1pt} {\kern 1pt}{\rm{ = }}{\kern 1pt} {\kern 1pt}{\gamma _2},}\\ {\tau {\kern 1pt} {\kern 1pt}{\rm{ = }}{\kern 1pt} {\kern 1pt}{\eta _{\rm{p}}}{\gamma _1} + {\eta _{\rm{s}}}{{\dot \gamma }_2} + \tau (H).} \end{array}} \right\}$

式中：η_p为虚数黏度，η_s为动态黏度，τ（H）为磁致剪切屈服强度.

2）串联黏弹塑性模型（串联模型）. 串联黏弹塑性模型由一个线性弹簧与一个Bingham黏塑性模型串联而成，如图3所示，本构关系为

(3) $\left. {\begin{array}{*{20}{l}} {\gamma {\kern 1pt} {\kern 1pt}{\rm{ = }}{\kern 1pt} {\kern 1pt}{\gamma _1}{\rm{ + }}{\gamma _2},}\\ {\tau {\kern 1pt} {\kern 1pt}{\rm{ = }}{\kern 1pt} {\kern 1pt}{\eta _{\rm{p}}}{\gamma _1} = {\eta _{\rm{s}}}{{\dot \gamma }_2} + \tau (H).} \end{array}} \right\}$

式中：γ₁为弹性应变，γ₂为黏性应变.

磁势F_M^[17]可以表达为

(4) ${F_{\rm{M}}} = NI = \phi {R_{\rm{m}}}.$

式中：R_m为磁阻，ϕ为磁通，N为线圈匝数，I为励磁电流.

磁通为 $\phi = \int\limits_S {B{\rm{d}}S}$，其中B为通过面积S的磁感应强度.

磁路的基本结构如图4所示. 图中，h为阻尼间隙间隙宽度，D为活塞直径，L为活塞长度，t为工作缸筒厚度，d为活塞杆直径，2b为活塞有 效长度.

对于图4所示的磁路结构，磁路形成闭合回路，且磁路各处磁通量应为恒值，根据磁路将活塞及缸体分为区域1、区域2、区域3，如图5所示. 阻尼间隙处的磁阻^[17]为

(5) ${R_{{\rm{MR}}}} = \frac{{\ln {\kern 1pt}\left( {1 + 2g/{D}} \right)}}{{2{\mu _{\rm{0}}}{\mu _{{\rm{MR}}}}{\text{π}}b.}}$

式中：μ _MR为磁流变液的磁导率.

对于区域1，磁阻^[17]为

(6) $ {R_{{\rm{m1}}}} = \frac{{L - 2b}}{{{\mu _0}{\mu _{{\rm{rp}}}} \cdot {{{\text{π}} }}\left[ {{{\left( {D - 2h} \right)}^2} - {d^2}} \right]/4}}. $

同理可知，区域2、3的磁阻^[17]分别近似为

(7) ${R_{{\rm{m}}2}} = \frac{{\left( {D - d} \right)/2}}{{{\mu _0}{\mu _{{\rm{rp}}}} \cdot {\text{π}}Db}},$

(8) ${R_{{\rm{m}}3}} = \frac{{L - b}}{{{\mu _0}{\mu _{{\rm{rh}}}} \cdot {\text{π}}\left( {D + g + t/2} \right)t}}.$

式中：μ_rp、μ_rh分别为活塞材料和工作缸材料的相对磁导率.

总磁势^[17]为

(9) ${F_{\rm{M}}} = \phi \left( {{R_{{\rm{m}}1}} + 2{R_{{\rm{m}}2}} + {R_{{\rm{m}}3}} + {R_{{\rm{MR}}}}} \right).$

根据式（4）、（9）可知，

(10) $ {{B_{{\rm{MR}}}} = \frac{{NI}}{{{R_{\rm{m}}}{S_{{\rm{MR}}}}}} = } {\frac{{NI}}{{\left( {{R_{\rm{m1}}} + 2{R_{\rm{m2}}} + {R_{\rm{m3}}} + {R_{\rm{MR}}}} \right){S_{{\rm{MR}}}}}}.} $

在动态剪切场下，磁流变液在屈服后期和屈服前期之间不断跳动. 由于磁流变液在屈服前具有黏弹性材料的特征，在正弦应变γ = γ₀sin (ωt)作用下，磁流变液的切应力和切应变之间的关系可以用复模关系表示：

(11) $ \tau = {G^*}\gamma . $

式中：G*为复切变模量，可以表示^[16]为

(12) $\left. {\begin{array}{*{20}{l}} {{G^*} = G' + G''i,}\\ {\tau = G'{\gamma _0}\sin {\kern 1pt} {\kern 1pt} (\omega t) + G''{\gamma _0}\cos {\kern 1pt} {\kern 1pt} {\kern 1pt} (\omega t).} \end{array}} \right\}$

式中：G'为储存模量，与单位体积材料变形储存能量的平均值成正比；G"为损失模量，与单位体积材料一个周期变形后消耗的能量成正比. 储存模量与损失模量的比称为损失因子^[16].

(13) $\delta (\omega ) = \frac{{G''(\omega )}}{{G'(\omega )}}.$

式中：ω为磁流变液的剪切加载角频率.

Ginderk等^[16]的研究表明，储存模量G'可以近似用下式计算：

(14) $G' \approx 3\varphi {\mu _0}{M_{\rm{s}}}H.$

式中：φ为磁性颗粒的体积分数，μ₀为真空磁导率，M_s为磁性颗粒的饱和磁化强度.

根据流变力学公式可知，在并联模型中磁流变液的附加刚度系数为k_p=pG'，其中p为与磁流变液工作状态有关的结构参数. 根据何小伟等^[18]推导出的并联等效刚度与串联等效刚度的关系可知，串联模型下的等效刚度为

(15) $ {k_{\rm{s}}} = \frac{{k_{\rm{s}}^{\rm{2}} + c_{\rm{s}}^{\rm{2}}{\omega ^2}}}{{c_{\rm{s}}^{\rm{2}}{\omega ^2}}}{k_{\rm{p}}}. $

磁流变阻尼器中存在弹性力、黏性力及摩擦力，在串联模型中阻尼器弹性力等于黏性力，因此在串联模型中阻尼器出力^[15]为

(16) $ F = {F_{\rm{f}}} + {F_{\rm{k}}} = {F_{\rm{f}}} + {F_{\rm{c}}}. $

式中：F为阻尼器总出力，F_k为弹性力，F_c为阻尼力，F_f为摩擦力.

根据平行平板间隙流动理论可知，阻尼器阻尼力^[16]为

(17) $\begin{split} {F_{\rm{c}}} = &\left( {\frac{{12\eta LA_{\rm{p}}^2}}{{\displaystyle{\rm{{\text{π}} }}D{h^3}}} + \frac{{{L\rm{{\text{π}} }}D\eta }}{h}} \right){u_{\rm{c}}}\left( t \right) + \\ &\left( {L{\text{π}} D\tau \left( {{B_{{\rm{MR}}}}} \right) + \frac{{3L\tau \left( {{B_{{\rm{MR}}}}} \right){A_{\rm{p}}}}}{h}} \right){\mathop{\rm sgn}} \left( {{u_{\rm{c}}}(t)} \right) = {F_\eta } + {F_\tau }. \end{split} $

式中：τ（B_MR）为磁流变液的磁致屈服强度，u_c（t）为磁流变液黏性流速，F_η为黏滞阻尼力，F_τ为库仑阻尼力.

阻尼器中的弹性力为

(18) $ {F_{\rm{k}}} = {k_{\rm{s}}}{x_{\rm{k}}}. $

式中：x_k为阻尼器工作过程中的弹性位移.

阻尼器可调范围为

(19) $ D = \frac{{{F_{{B_{\max }}}} - {F_{{B_0}}}}}{{{F_{{B_0}}}}} = \frac{{{F_{\tau \left( {{B_{\max }}} \right)}}}}{{{F_\eta } + {F_{\rm{f}} }}}. $

式中： ${F_{{B_{\max }}}}$为阻尼器磁饱和时的最大出力， ${F_{{B_0}}}$为阻尼器的零场阻尼力.

如图6所示为图7所示剪切阀式磁流变阻尼器的结构示意图. 该阻尼器由活塞杆、活塞、磁流变液、缸体、端盖、密封环及线圈组成. 励磁单元由绕制在活塞上的线圈形成的电磁体组成，磁力线经由活塞、磁流变液和缸体构成回路.

在工作状态下，活塞杆带动活塞挤压腔体内的磁流变液，上、下腔体之间出现压力差，使得磁流变液通过阻尼间隙由高压区流向低压区. 当磁流变液流过阻尼间隙时，通过改变电流调节垂直于流体流动方向的磁感应强度，使得磁流变液发生流变效应，实现阻尼力连续控制.

因为该剪切阀式磁流变阻尼器阻尼器为回转类装置，为了简化模型，可以采用二维轴对称模型进行磁流变阻尼器的磁场仿真分析. 该阻尼器的活塞、磁流变液等材料具有非线性的磁化特性，所以需要在COMSOL中导入图8所示各材料的B-H曲线. 在阻尼器工作过程中，活塞带动线圈在缸体内不停地上下运动，所以在该过程中的磁场是不断变化的. 为了仿真运动过程中的阻尼器磁场分布，增加了变形几何，给活塞与线圈设定指定的位移，通过仿真计算阻尼器运动状态下的磁场.

通过变形几何，给活塞杆、活塞和线圈施加一个频率为2 Hz、振幅为3 mm的正弦位移，给线圈通1.2 A的电流，t=0 s时的磁场仿真结果如图9所示.

t=0.125 s时的磁场仿真结果如图10所示.

根据图9、10可知，磁感应强度分布随着活塞与线圈的运动而变化，且阻尼通道内的磁感线一直集中在活塞的两侧.

根据式（2）可知，并联模型下磁流变液的表观黏度为

(20) $ \eta \,=\, \frac{\tau }{{\dot \gamma }}\,{\rm{ = }}\,{\eta _{\rm{p}}}\frac{\gamma }{{\dot \gamma }} + \frac{{\tau (H)}}{{\dot \gamma }} + {\eta _{\rm{s}}}. $

根据式（3）可知，串联模型下磁流变液的表观黏度为

(21) $ \left. {\begin{array}{*{20}{l}} {\eta \,{\rm{ = }}\,\tau /\dot \gamma ,}\\ {\tau = {\eta _{\rm{s}}}\dot \gamma + \tau (H) - \dfrac{{{\eta _{\rm{s}}}}}{{{\eta _{\rm{p}}}}}\dfrac{{\partial \tau }}{{\partial t}}.} \end{array}} \right\} $

根据图9、10可知，当线圈通电时，阻尼通道中的磁场主要分布在活塞附近的磁流变液上. 根据图11磁场与磁流变液剪切屈服强度的关系可知，磁场越大，磁流变液的剪切屈服强度越大，所以活塞附近的磁流变液剪切屈服强度远远大于其他区域. 将该区域称为激活区，其他区域称为未激活区，如图12所示.

将磁感应强度分布数据导入三维流固耦合场，如图13所示. 根据Lord公司提供的材料参数，得到磁流变液的磁致剪切屈服强度分布，如图14所示. 通过偏微分方程建立磁流变液的应力-应变关系，利用幂律模型设置流体表观黏度. 将缸体与上、下法兰固定，给活塞杆和活塞施加正弦位移激励，使得上、下腔体内的磁流变液通过阻尼通道由高压区流向低压区，通过仿真得到在磁流变液流动时阻碍活塞运动的反作用力.

如图15所示，使用直流电源给阻尼器线圈通电，调节阻尼器的磁场并改变阻尼器出力，通过液压执行器带动磁流变阻尼器活塞杆作简谐运动，通过数据采集系统采集力-位移信号. 根据不同的磁流变液本构关系对磁流变阻尼器的出力进行仿真，将仿真数据与示功实验数据进行对比，结果如图16~19所示. 图中，f为频率，A为振幅，I为电流.

对比图16~19的仿真与实验示功曲线可以看出，在不同的f、A和I下，根据串联本构建立的有限元仿真模型得到的阻尼器出力结果与示功实验结果吻合度更高. 综合图16~19可以看出，在阻尼器工作过程中电流对阻尼器出力的影响最大. 根据图18、19可以看出，改变电流，阻尼器的刚度会随之改变，且变化明显.

在频率为4 Hz、振幅为0.5 mm的正弦激励下，通过理论仿真分析串联模型下不同结构参数对阻尼器在电流为0 A和2 A时的最大出力F_max及其可调范围的影响.

图20中，∆α为结构参数的变化比例.根据图20可知，不加电流时阻尼器出力随着活塞杆直径和阻尼间隙的增大而减小，随着活塞直径和活塞有效长度的增大而增大. 根据图21、22可知，阻尼器出力与可调范围都随着活塞杆直径的增大而减小，随着活塞直径的增大而增大，随着活塞有效长度与阻尼间隙的增大阻尼器出力都有先增大后减小的趋势. 根据图20~22可知，活塞直径对阻尼器出力与可调范围的影响最大. 当活塞直径增大50%时，阻尼器最大出力为1 496 N，阻尼力可调倍数为9. 33.

灵敏度函数描述参数对阻尼器出力影响随时间的动态变化，不能量化参数变化对阻尼力变化的影响程度，因此采用峰值灵敏度指标与均值灵敏度指标^[19]量化阻尼器各结构参数对阻尼器出力的影响.

峰值灵敏度指标为

(22) ${S_1} = \frac{{\left| {\Delta {F_i}} \right|}}{{{F_{\max }}}}\left| {\begin{array}{*{20}{c}} {} \\ {_{\max }} \end{array}} \right. \times 100{\text{%}} = \frac{{\left| {{\lambda _{\alpha i}}} \right| \cdot \Delta {\alpha _i}}}{{{F_{\max }}}}\left| {\begin{array}{*{20}{c}} {} \\ {_{\max }} \end{array}} \right. \times 100{\text{%}} .$

均值灵敏度指标为

(23) ${S_2} = \int_0^{{t_0}} {\left| {{\lambda _{\alpha i}}} \right| \cdot \Delta {\alpha _i}{\rm{d}}t} .$

如图23、24所示为当结构参数减小50%时，各结构参数的峰值与均值灵敏度指标. 根据峰值灵敏度指标与均值灵敏度指标计算可知，活塞直径与阻尼间隙对阻尼出力的影响较大，其中活塞直径的影响大于阻尼间隙，当结构参数减小50%时活塞直径的峰值灵敏度指标与均值灵敏度指标分别为84.66%和94.51 N.

（1）根据有限元仿真建立基于复杂流体本构的阻尼器出力模型，通过实验验证该模型的准确性，本文提供了复杂本构流体的模拟方法. 该方法能够较准确地模拟磁流变液复杂的流变效应，在设计阶段较准确地模拟阻尼器出力性能，为阻尼器的结构设计提供依据.

（2）通过理论与仿真比较发现，相较于其他3个结构参数，活塞直径对阻尼器最大出力及可调倍数的影响最大，当活塞直径增大50%时阻尼器最大出力为1 496 N，阻尼力可调倍数为9.33.

（3）通过灵敏度分析可知，阻尼间隙与活塞直径对阻尼器出力的影响较大，主要是由于在剪切阀式磁流变阻尼器工作过程中，磁流变液压差流动引起的反作用力大于剪切流动引起的反作用力.