马良^1,2, 肖萍², 王任胜^1,3, 修世超¹
1. 东北大学 机械工程与自动化学院, 辽宁 沈阳 110819;
2. 新疆工程学院 机械工程系, 新疆 乌鲁木齐 830091;
3. 辽宁科技学院 机械工程学院, 辽宁 本溪 117004
收稿日期：2016-09-12
基金项目：国家自然科学基金资助项目(51375083);沈阳市科学技术计划项目(F16-205-1-02);辽宁省教育厅科学研究项目一般项目(L2015274)。
作者简介：马良(1976-),男,甘肃临洮人,东北大学博士研究生;
修世超(1958-), 男, 辽宁凌源人, 东北大学教授, 博士生导师。

摘要：为降低采用标准孤立链理论模型计算高体积分数磁流变液流变剪应力时的计算误差, 提出一种通过距离权重系数评价磁链间相互作用的流变模型.编写Monte Carlo仿真程序计算颗粒位置与磁矩方向矢量、同时输出与系统颗粒间距有关的距离权重系数.按标准和距离权重这两种链结构模型计算磁流变液剪切应力数值, 并和4种磁流变液(30.58 % ~44.79 %)的流变仪剪应力测试数据进行比较, 结果显示, 权重模型在磁场全范围内具有更高的计算精度, 并且和流变试验实际结果基本符合, 证实了磁链微观结构间的相互作用力是影响磁流变液流变特性的主要因素之一.
关键词：磁流变液流变特性距离权重系数剪切屈服应力孤立链模型蒙特卡洛方法
An Isolated Chain Rheological Shear Stress Model with Distance Weighted Coefficients
MA Liang^1,2, XIAO Ping², WANG Ren-sheng^1,3, XIU Shi-chao¹
1. School of Mechanical Engineering & Automation, Northeastern University, Shenyang 110819, China;
2. Department of Mechanical Engineering, Xinjiang Institute of Engineering, Urumqi 830091, China;
3. School of Mechanical Engineering, Liaoning Institute of Science and Technology, Benxi 117004, China
Corresponding author: MA Liang, E-mail:maliang7653@sina.com
Abstract: For reducing the error in calculating shear stress of high volume concentration magnetorheological fluids by applying standard isolated chain theoretical model, a new rheological model was proposed to evaluate interaction of magnetic chains by distance weighted coefficients.A Monte Carlo program was written for simulating position coordinates and direction of magnetic moment of carbonyl iron particles in magnetorheological fluids system, with two distance-related coefficients output.Comparing theoretical values of standard and distance weighted shear stress model with four magnetorheological fluids' testing data that measured by rheometer, results indicated that distance model has higher computational accuracy at the whole experimental range of magnetic field than that using standard isolated chain-based one, and it generally aligned with rheological results, which proves interaction force between magnetic chains is one of impact factors that can significantly affect rheological characteristic of magnetorheological fluids.
Key Words: magnetorheological fluidrheological characteristicdistance weighted coefficientshear yield stressisolated chain modelMonte Carlo method
磁流变液(magnetorheological fluids, MRFs)是新型智能材料, 由载液、分散悬浮相铁磁颗粒(carbonyl iron particles, CIPs)、防止沉降的稳定剂等构成.对MRFs施加外磁场, 颗粒沿磁场方向聚集排布, 固相抗剪强度比零磁场条件下高出2~5个数量级^[1].这一性质广泛应用于磁流变液联轴器、磁流变阀、离合器、磁流变抛光等机械设备与加工技术, 因此近年来, 发生流变效应后对MRFs磁致剪应力的计算分析, 是材料和机械加工领域的热点问题之一.
一些微观剪力模型把场致微结构假设为链、层或板形结构, 考虑颗粒极化作用, 建立平衡方程得到磁致特性参数^[2-3].物理试验也证实了低体积分数特定条件下模型结果的准确性^[4-5].但颗粒成链受外磁场、空间位阻以及热运动等因素影响, 成链后颗粒坐标信息无法确知, 如果假定系统处于某种单一形式的构造体系, 当颗粒体积分数较高时, 存在较大的误差^[6].
为缩小结构假设误差, 分析过程分成两步：首先建立孤立链剪力微结构计算模型, 对磁流变液单元系统进行Monte Carlo仿真, 输出施加外磁场后的颗粒坐标和磁矩方向, 并计算两两颗粒的间距；再按间距判断所有颗粒的成链条件, 建立有效磁链集合, 遍历集合内每条磁链, 基于前一步骤的仿真位形信息得到距离权重剪应力模型的计算值.
1 孤立链剪切应力理论模型如图 1所示：剪切力τ的作用区域是底面积为D²、高h=λ_iR的六面体, 其中D=λ_sR；λ_i, λ_s分别为颗粒中面垂直和平行于B₀的径向距离系数.如果颗粒体积分数为φ, 则λ_s, λ_i, φ三者间按作用区域的体积等价关系满足：

(1)

图 1(Fig. 1)

图 1 剪力计算受力分析简图Fig.1 Shear force analysis diagram

如图 1, 由安培分子电流假说和偶极子理论, 磁化作用力F是以O₀上磁化电流对任意点P的作用为被积表达式, 在颗粒O_i表面的累积和：

(2)

式(2)颗粒感应强度被积表达式根据偶极子理论写成式(3)的形式.其中：μ_f为基液磁导率; m₀为磁偶极矩, r为磁偶极子中心指向P点的矢量.

(3)

颗粒在水平力作用下从虚线位置移动至O_i位置, 颗粒间剪切角为ω, 中心距扩大为h^ω=h/cosω.如图 1, 式(2)作用力表达式通过坐标变换得到沿x轴和z轴方向的两个分力F_ix和F_iz, 通过受力分解, 水平方向作用力的合力即所求F_κ：

(4)

F_ix和F_iz由式(2)中颗粒O₀对O_i上点P处磁感应强度B_P积分式决定, 联立式(1)~式(4), 经坐标变换^[7], 简化得

(5)

系数κ(λ_i, ω)按式(6)计算, 其中f_x(λ_i, ω)和f_z(λ_i, ω)是F_ix和F_iz的两个系数, 包含三重积分且没有解析解, 用MATLAB数值积分命令integral3计算.

(6)

选择Fr?elich-Kennelly材料模型^[8], 颗粒磁化强度如式(7)所示.其中χ₀为铁磁颗粒材料的磁化率, 取χ₀=1 000；H为外磁场强度(kA/m)；M_s为颗粒磁化饱和强度.

(7)

式(5)~式(7)就是与距离系数λ_i和颗粒磁化强度M有关的孤立链微观剪应力计算表达式.
式(5)中磁链沿水平面方向的控制区域D²和颗粒体积分数φ有关, 通常均匀扩散的磁流变液外覆稳定剂层的厚度中位数δ为半径的10 %, 故设相邻颗粒成链极限间距为1.1d=2.2R.
令纵向距离系数λ_i=λ_s=2.2, 代入式(1)解得磁链相互不影响的极限值为φ=39.34 %, 此时外磁场水平或垂直两个方向均满足成链条件, 如图 2所示.颗粒磁矩m和外场B₀存在夹角, 即水平方向分量不为零时, 磁链间有聚集趋势, 颗粒间作用力不再只沿外场方向, 而是空间相互叠加作用的综合结果, 因此按单链模型计算, 理论结果比实际数值偏小.
图 2(Fig. 2)

图 2 磁链在垂直外磁场方向的相互影响Fig.2 Chains' interaction effects in the direction vertical to exterior magnetic field

为降低孤立链模型在颗粒体积分数大时的计算误差, 把磁链横向间距的缩小转为用颗粒磁化强度的增幅进行表述.通过Monte Carlo仿真位形, 从有效磁链中提取位形信息并计算颗粒距离和等效磁化强度增幅因子, 以权重值代替孤立链模型表达式(5)中的λ_i和M, 计算剪应力修正值.
1.1 Monte Carlo仿真输出位形矩阵定义边长的立方体单元系统：N个直径d的球形颗粒在初始位形均布于立方体空间, φ是铁磁颗粒体积分数.
对微米级粒径磁流变液, 忽略热运动对系统的影响, 故磁性颗粒能量主要由偶极子极化能E_ij^dd、空间位阻能E_ij^s和外磁场在颗粒位置的势能构成：

(8)

为避免当颗粒间距r_ij =| r_i-r_j| < r_ij^min时, 位置相互“侵入”的情况, 设定此时空间位阻能为无穷大, 因此前后两次循环总势能计算值之差ΔE⁽ⁿ⁾?0, 这相当于一个惩罚性高“势垒”, 以此来拒绝错误位形.
成链颗粒中心距r_ij^min既要保证颗粒外覆稳定剂存在, 又不至间距过大导致极化能急速衰减而无法成链.综合考虑之后, 磁链上颗粒“相邻”的判定依据是颗粒中心距在粒径1.1~1.2倍范围内.
颗粒i因磁化对颗粒j产生的极化能E_ij^dd如式(9)所示, 其中:r_ij为颗粒中心距矢量；μ₀为真空磁导率, μ₀=4π×10^-7 N/A²; m_i为颗粒i的磁矩, ^[9], M_s为颗粒磁化饱和强度.

(9)

Monte Carlo仿真构造足量颗粒的随机移动样本集, 按照式(8)和式(9), 结合“总势能下降”和“几何坐标不冲突”两个约束条件, 循环判断是否接受颗粒移动.直至相邻两次循环能量差ΔE_i⁽ⁿ⁾减小到稳定数值, 输出末状态位形坐标矩阵P₀：

(10)

1.2 有效磁链判断按中心距r_ij满足r_ij∈[1.1d, 1.2d]的条件, 任选P₀集合中索引编号为1的颗粒, 寻找与其满足距离要求的颗粒i, 如有多个颗粒同时满足距离要求, 如图 3右侧下方r₁ < r₂的两个颗粒, 取其中距离最小的i替换原颗粒1, 重复上述步骤直至没有颗粒与之满足相邻, 判定构成第i条颗粒数为k_i的磁链, 记为C_i.
图 3(Fig. 3)

图 3 颗粒的磁链归属性判断Fig.3 Judgement of chain belonging of CIPs

对颗粒数N=64、体积分数38.7 %的单元系统施加H=250 kA/m的外磁场, Monte Carlo仿真末态位形微结构存在如图 4所示的四类基本链结构.分析发现：磁链(图 4b和图 4c)受剪切作用拉伸并吸附颗粒时会重新构成如图 4a所示的典型有效直链；而磁链(图 4d)尽管多数颗粒磁矩方向与外场方向一致, 但相邻颗粒中心距与外场过大的夹角决定了相邻颗粒间吸引作用微弱, 分离趋势更大, 其上颗粒更容易吸附于其他磁链并重组, 不能作为有效磁链.
图 4(Fig. 4)

图 4 系统中不同的磁链类型Fig.4 Different kinds of chains in MRFs system (a)—直链(有效链); (b)—L形链(有效链)； (c)—含吸附颗粒的直链(有效链); (d)—直链(无效链).

综上, 编程计算各磁链相邻颗粒中心距矢量r_ij与外场间夹角均值, 用夹角作为有效磁链判断的阈值^[10], 如夹角均值大于该角度则判为无效磁链.最终形成有效磁链集合{C_i} _{i=1, …, n}
1.3 按位形权重修整计算模型参数单元子系统内n条有效磁链的距离权重系数均值记为{λ_Ci}|_{i=1, 2, …, n}, 第C_i条磁链上, 颗粒总数记为k_i, 则该磁链上距离权重系数均值为

(11)

式中λ_Ci(t)为磁链C_i上第t个颗粒对的距离权重系数.同时提出权重增幅系数A_M, 用于评价磁链之间的相互作用, A_M与系统颗粒总数N、颗粒块体材料属性调整系数, 以及系统内所有n条有效磁链上颗粒数量的均值k有关.

(12)

用式(11)中的λ_Ci、式(12)得到的增幅磁场强度A_M·M代替式(5)中的λ_i和M, 得到含有距离权重的流变剪应力计算模型：

(13)

2 计算结果与实验验证自行制备水基磁流变液, 用十二烷基苯磺酸钠(体积分数2.2 %)加入去离子水(体积分数59.1 %), 加热至溶解, 放入饱和磁化强度M_s=1.4T的JCF1-2型羰基铁粉(体积分数38.7 %, 吉林吉恩镍业公司)搅拌1 h至无沉淀^[11].图 5为磁流变液试样中外覆稳定剂铁磁颗粒的电镜扫描图, 可以看出粒径分布约为2~3 μm, 同时出现明显絮状稳定剂层, 综合权衡后, 取测试样品平均粒径为d=2.2 μm.
图 5(Fig. 5)

图 5 磁流变液试样中铁磁颗粒SEM图片Fig.5 CIPs' SEM image in MRFs sample

将制备的磁流变液试样略干燥, 置于图 6所示流变仪固定板, 测定0~250 kA/m磁场强度变化范围内, 不同剪切速率下试样的剪应力τ值.
图 6(Fig. 6)

图 6 MCR301流变仪及测量原理Fig.6 MCR301 rheometer and measuring principle

从图 7看出, 不同磁场强度下样品的剪应力变化趋势基本相同, 且剪切速率在[0.01, 1]范围是剪力攀升阶段, 剪切速率达到[1, 100]区间时, 剪力基本维持恒定, 可以认为达到了该磁场强度条件下的剪应力稳定值.
图 7(Fig. 7)

图 7 不同磁场强度下的剪切应力实验值Fig.7 Experimental values of shear stress at different magnetic field intensities

取剪切速率为10 s^-1时的剪应力测试数据τ_exp, 调用Monte Carlo程序得到距离权重系数, 以式(13)作为磁感应强度B的自定义函数τ_B, 将调整系数β的计算视作对数据组(B, τ_exp)寻找最佳拟合的最小二乘拟合问题, 通过MATLAB中的lsqnonlin指令计算, 代入式(12)得到A_M, 上述数据列入表 1.
表 1(Table 1)

表 1 颗粒磁化强度权重系数Table 1 Weighted factor of CIPs' magnetization intensity H/(kA ·m-1) 82 142 186 220 250 k 5.31 4.23 4.25 5.44 5.11 λCi 2.241 2.245 2.206 2.246 2.232 β 0.250 6 0.232 9 0.201 2 0.258 4 0.239 6 AM 1.273 1 1.277 4 1.235 1 1.279 6 1.264 0

表 1 颗粒磁化强度权重系数 Table 1 Weighted factor of CIPs' magnetization intensity

为进一步验证权重剪应力模型对不同材料和体积分数磁流变液的适应性, 增加三种LORD公司制备的合成油基磁流变液的τ-B实验数据, 型号依次为MRF-122EG, MRF-132DG和MRF-140CG, 物理特性列于表 2.
表 2(Table 2)

表 2 磁流变液物理特性参数Table 2 Typical properties parameters of MR fluids 型号 ρ/(g·cm-3) 体积分数/% M1/T MRF-122EG 2.28~2.48 30.58 1.171 8 MRF-132DG 2.95~3.15 39.57 1.407 1 MRF-140CG 3.54~3.74 44.79 1.565 7

表 2 磁流变液物理特性参数 Table 2 Typical properties parameters of MR fluids

用B=μ(H+M)模拟颗粒非线性磁化过程, 设磁流变液整体饱和磁化强度M_s与颗粒块体饱和磁化强度M间满足M_s=φM.M数据按产品页面公布的H-B曲线数据, 以Fr?elich-Kennelly模型做最小二乘拟合, 图 8为MRF-122EG模型MRFs拟合结果与按照B=μ₀H的线性磁化计算结果以及H-B实验数据的比较.
图 8(Fig. 8)

图 8 MRF-122EG磁化曲线(M=1.1718 T)Fig.8 Magnetization process curves of MRF-122EG

从图 8看出, Fr?elich-Kennelly模型对颗粒材料磁化过程的描述与实验数据基本吻合.
最后, 按标准单链模型、权重模型和实验数据绘制4种磁流变液的屈服剪切应力τ_y与磁感应强度B的关系曲线, 见图 9.观察发现：对于体积分数φ∈[30.58 %, 44.79 %], 且不同载液类型的4种磁流变液, 应用两种剪应力微观模型得到的应力变化趋势相同, 随颗粒磁感应强度增加而尚未达到磁化饱和强度前, 孤立链模型计算剪应力增长速率远低于实验数据, 权重模型计算值则在颗粒磁感应强度几乎全范围内与实验数据较为符合.
图 9(Fig. 9)

图 9 模型计算值和实验结果比较Fig.9 Comparison of calculating values of two models and experimental data (a)—MRF-122EG, φ=30.58 %; (b)—MRF-132DG, φ=39.57 %； (c)—MRF-140CG, φ=44.79 %; (d)—自备样品, φ=38.70 %.

模型计算值对比实验值, 拟合优度R²如表 3所示.可看出孤立链模型R²值均低于0.5, 说明孤立链模型不适合描述高体积分数MRFs的流变特性；比较结果证实微观形态磁链间的确存在相互作用, 且对磁流变液剪应力变化有显著影响.
表 3(Table 3)

表 3 剪应力计算值的R2Table 3 R2 of shear stress calculating data 型号 孤立链模型 距离权重模型 水基制备MRFs 0.471 9 0.969 9 MRF-122EG -1.354 7 0.943 2 MRF-132DG -1.850 7 0.897 8 MRF-140CG -2.671 2 0.860 3

表 3 剪应力计算值的R2 Table 3 R2 of shear stress calculating data

3 结语高体积分数磁流变液普遍用于机械领域的磁流变抛光液制备, 但当颗粒体积分数增大时, 孤立链剪应力微观计算模型因未考虑磁场控制区域重叠造成计算误差显著增大, 就此问题提出一种结合Monte Carlo仿真位形的修正计算方法.这种方法弥补了孤立链计算模型一些与几何位形相关的参数无法测定的问题, 更重要的是, 在孤立链模型中引入磁链颗粒权重经验增幅系数, 用于描述聚集磁链间因磁性颗粒局部磁场作用区域叠加所导致的磁场增幅效果, 使得权重模型适合于解释不同载液环境下, 较高体积分数的磁流变液流变剪切应力的变化规律.通过拟合优度R²对比不同磁场强度下, 孤立链模型计算值和流变仪测试结果, 证实了如下结论：权重模型显著提高了对高体积分数磁流变液的流变剪切应力预测精度, 且计算模型形式简单, 适合在实际工程中对磁流变液的流变力学特性进行预测与分析.
参考文献

[1]	Bossis G, Iskakova L, Kostenko V, et al. Kinetics aggregation of magnetic suspensions[J].Physica A:Statistical Mechanics and Its Applications, 2011, 390(14): 2655–2663.DOI:10.1016/j.physa.2011.02.044
[2]	Weis J J, Levesque D. Chain formation in low density dipolar hard spheres:a Monte Carlo study[J].Physical Review Letters, 1993, 71(17): 2729–2732.DOI:10.1103/PhysRevLett.71.2729
[3]	Tang X, Zhang X, Tao R, et al. Structure-enhanced yield stress of magnetorheological fluids[J].Journal of Applied Physics, 2000, 87(5): 2634–2638.DOI:10.1063/1.372229
[4]	Kruse T, Krauth?user H G, Spanoudaki A, et al. Agglomeration and chain formation in ferrofluids:two-dimensional X-ray scattering[J].Physical Review:B, 2003, 67(9): 094206.DOI:10.1103/PhysRevB.67.094206
[5]	Rodríguez-López J, Castro P, Elvira L, et al. Study of the effect of particle volume fraction on the microstructure of magnetorheological fluids using ultrasound:transition between the strong-link to the weak-link regimes[J].Ultrasonics, 2015, 61(5): 10–14.
[6]	高伟, 魏齐龙, 李晓媛, 等. 磁流变抛光与磁流变液:原理与研究现状[J].磁性材料及器件, 2015, 46(2): 68–73. ( Gao Wei, Wei Qi-long, Li Xiao-yuan, et al. Magnetorheological finishing and magnetorheological fluid:principle and research situation[J].Journal of Magnetic Materials and Devices, 2015, 46(2): 68–73.)
[7]	Zhu C C. A new theoretical model about shear stress in magnetorheological fluids with small shear deformation[J].Journal of Wuhan University of Technology(Materials Science), 2005, 20(1): 52–56.DOI:10.1007/BF02870873
[8]	Ivaneyko D, Toshchevikov V, Saphiannikova M, et al. Effects of particle distribution on mechanical properties of magneto-sensitive elastomers in a homogeneous magnetic field[J].Condensed Matter Physics, 2012, 15(3): 33601–33612.DOI:10.5488/CMP.15.33601
[9]	Bradbury A, Menear K, O'Grady K, et al. Magnetic size determination for interacting fine particle systems[J].IEEE Transactions on Magnetics, 1984, 20(5): 1846–1848.DOI:10.1109/TMAG.1984.1063338
[10]	李海涛, 彭向和, 黄尚廉. 基于偶极子理论的磁流变液链化机理模拟研究[J].功能材料, 2008(6): 902–904. ( Li Hai-tao, Peng Xiang-he, Huang Shang-lian. Study on the chain-formation mechanism of magnetorheological fluids based on dipole theory[J].Journal of Functional Materials, 2008(6): 902–904.)
[11]	孙百万, 黎胜权, 王任胜, 等. 往复式动磁场磁流变抛光机理及抛光液制备[J].机械设计与制造, 2016(7): 81–84. ( Sun Bai-wan, Li Sheng-quan, Wang Ren-sheng, et al. Mechanism of the magnetic rheological polishing of reciprocating dynamic magnetic field and preparation of polishing liquid[J].Machinery Design & Manufacture, 2016(7): 81–84.)