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

摘要：热运动对纳米级磁流变液系统的颗粒团聚行为和宏观流变特性扰动影响不可忽视.为研究施加外磁场后系统微结构与热运动间的关系, 探讨在不同颗粒粒径情况下热运动能量在系统总能量中的占比.采用CCD设计方法, 分析Monte Carlo仿真获取的颗粒位形信息, 生成关于外磁场磁感应强度、颗粒粒径和体积分数的系统热耦合系数 <λ>回归模型; 成链团聚的临界颗粒粒径计算值比文献经验值高出约23%, 证实在考虑系统内多因素的综合影响下, 纳米级磁流变液具有更严格的流变现象临界发生条件.
关键词：磁流变液热耦合系数中心组合设计试验设计流变特性
Chain Formation and Agglomeration Characterization
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: Thermal motion behavior has an unneglectable disturbing effect on the agglomeration and rheological characteristics of nano-size magnetorheological fluids(MRFs).To discover the relationship between microstructure and thermal motion in MRFs with external magnetic fields applied, the thermal proportion of total energy with different diameters of carbonyl iron particles (CIPs) was discussed. The CCD method was employed to analyze the CIPs' position data set calculated by Monte Carlo simulation, and a systematic thermal coupling coefficient <λ> regression model related to magnetic induction intensity, CIPs' diameter and volume concentration was built. The threshold particle diameter value calculated by the model is 23% higher than the empirical value in the literature, which proves that rheological characteristic has a stricter precondition with multiple factors considered in MRFs.
Key Words: magnetorheological fluids(MRFs)thermal coupling parametercentral composite design(CCD)experiment designrheological characteristic
磁流变液(magnetorheological fluids, MRFs)是一种优秀的智能材料, 其剪切模量、屈服强度等力学特性并非固定不变, 外磁场下毫秒量级时间内, 能从低黏度牛顿流体迅速转变为半固态低流动性、高黏度Bingham塑形体^[1-2], 控制磁化电流参数, 能获得远高出零磁场时的可调剪应力, 即“流变”效应.磁流变液的流变特性与微观结构参数间的关系引起国内外学者几十年来持续关注, 其中磁性流体在何种条件下产生流变效应, 是其中关键问题之一.
流变效应与磁性颗粒的团聚成链密切相关, 因此当磁性颗粒分布粒径达到纳米尺度时, 热运动程度将是场致磁性颗粒成链的潜在扰动因素之一.以往研究中, 李海涛等^[3]计算出常温条件下, 磁流变液(d=5 μm)内单独一对颗粒间的磁热能比值约为5.35×10⁶, 因此对微米量级磁流变液, 热运动能对成链特性的影响可忽略不计; 但纳米级磁流变液中情况有很大不同, 颗粒粒径的急剧降低导致热作用效应对颗粒运动干扰明显.此外, 从目前研究成果看, 粒径分布和颗粒团聚二者的具体关系还存在一定分歧, 例如Kormann等^[4]通过流变试验证实平均粒径30 nm的磁流变液动态屈服强度仍能达到τ=2.5 kPa (φ=23%, B=0.2 T); 小粒径颗粒布朗力对流变特性有较强扰动的研究时有报道^[5-6].因此工程制备磁流变液, 常使用高于100 nm的磁性颗粒, 以保持较大的磁热能比, 避免无法成链的情况^[7].但磁流变液属于固-液两相材质, 影响流变特性的因素有多种, 且诸因素可能存在交互影响, 目前尚未发现试验观测获取发生流变效应时临界参数取值的相关报道.
为准确评估MRFs系统整体的团聚成链, 以及流变特性发生的临界条件, 通过响应面分析和中心组合设计等统计方法, 以热耦合系数为响应值, 探讨颗粒体积分数、粒径和磁场强度3个主项因子的取值, 包括参数间存在的交互影响项对颗粒团聚贡献的量化描述方法.
1 MRFs系统颗粒团聚的判定准则1.1 热耦合系数引入无因次热耦合系数λ, 即:极化能与热运动能之比评价颗粒对之间的热运动幅度, 见式(1)^[8]:

(1)

其中:Boltzmann常数k_B=1.38×10^-23J·K^-1; T为温度; E_ij^dd是两颗粒间极化能, 代表位置r_i处磁矩方向矢量m_i的颗粒i, 受到r_j处磁感应强度为B_j(r_i)的颗粒j的磁偶极能作用, 如式(2)所示:

(2)

式中:磁矩模表示为|m_i|=M_d·πd³/6, M_d为颗粒平均磁化饱和强度, 磁矩单位向量定义为m_i/|m_i|; 真空磁导率取μ₀=4π×10^-7H/m; r_ij=r_j-r_i为两颗粒间的距离矢量; r_ij=|r_j-r_i|为颗粒间距.
从式(2)知E_ij^dd∝d³, 当粒径降低, 颗粒间极化束缚力随之以立方幂律降低, 热运动能占比相应变大, 颗粒越容易摆脱磁链吸附, 成链团聚越困难, 形成各向异性链式排布的概率越低, 流变剪切应力减小, 越不容易发生流变效应.此外, 宏观剪应力大小与颗粒间距直接相关, 但受外磁场分布、颗粒热运动、稳定剂外覆层阻隔、颗粒体积分数及粒径分布等多种因素影响, 磁致微结构的各向异性链式聚集形态复杂, 尚无文献表明实际试验可取得颗粒实时位形坐标, 也缺乏对磁流变液系统整体热耦合系数和流变力学参数间关系评估的模型.因此提出采用Monte Carlo仿真获得系统势能下降趋于稳态时的颗粒坐标与磁矩方向矩阵, 获取颗粒间距和磁矩方向信息^[9], 进而判断颗粒团是否聚成链的计算思路.
1.2 颗粒成链极限距离的约定利用热耦合系数判定相邻颗粒的团聚, 需要知道颗粒处于“聚集”或“成链”的极限间距.显然:理想条件下, 对2个等径颗粒, 该极限间距值应为r_ij^min=R_i+R_j=2R, 但实际上为防止颗粒沉降, 磁流变系统内将掺入一定比例的稳定剂防止颗粒沉降.考虑外磁场下颗粒的外覆稳定剂层以及层间部分渗透, 等径颗粒系统最小间距取为r_ij^min=1.1d; 此外, 根据式(2)可知, 间距变大时, 颗粒极化能急剧衰减, 给定颗粒成链的最大相邻间距为r_ij^chain=1.2 d^[10].于是颗粒“团聚”或处于“成链”的距离范围是

(3)

1.3 颗粒间的范德华势能颗粒间范德华能常用于描述磁流变液颗粒的热作用效果, 见式(4).

(4)

其中, A_H为Hamaker常数, 与颗粒介电属性、外覆稳定剂层和载液有关, 通常取A_H=0.5×10^-19J ^[11].比较式(2)和式(4)发现:颗粒间极化能和范德华能都与间距r_ij大小呈正比例相关, 但前者有远高于后者的增幅, 意味着范德华能对系统能量的影响与颗粒粒径有关; 对微米级磁性颗粒, |E_ij^vdw|＜＜|E_ij^dd|, 范德华能对系统位形和能量组成的影响可忽略不计; 当λ→1, 即:热运动能量和偶极子能间接近平衡, 此时分布粒径早已达到纳米量级, 范德华能对磁性颗粒能量组成的影响不能忽视.
1.4 Monte Carlo算法与流程对由N个颗粒组成的MRFs系统, 每次循环随机选择其中1个颗粒, 沿随机磁矩方向移动小的步长, 按式(5)所示空间位阻能定义, 以无穷大势垒阻止颗粒间不符合物理实际的相互“侵入”, 排除与系统其他颗粒的几何位置冲突^[12].

(5)

联立式(2), 式(4)和式(5), 按式(6)计算系统总位形能, 与前次循环的计算结果相减得位形能改变量ΔE_i.

(6)

根据Metropolis算法, 随机试探移动是“任意”的, 但是否接受本次试探移动, 则基于系统能量状态迁移概率判断, 如式(7)所示:

(7)

第i次循环玻尔兹曼因子是位移沿总势能下降的接受概率, ΔE>0时, 与本次试探位移r_ij⁽ⁿ⁾接受概率P_ij呈反比变化关系.颗粒移动时, 平动和旋转步长因子应保证位形接受概率约为30%~50%, 基本覆盖全部颗粒样本^[13].第n次循环时, 一系列接受位形概率P_ij构成空间Markov链, 相邻两次位形能差ΔE_i逐渐减小并收敛到某一数值, 输出能量组成及系统末状态位形坐标与磁矩方向.
由能量条件控制颗粒的接受位形应保证非常高的样本容量, 因此系统位形初值具体是何种形状对最终位形状态的影响可忽略不计.为简化计算起见, 初始条件考虑采用颗粒坐标矢量r_i按行列均匀规则排布的立方体结构, 边长L=.其中:N为系统颗粒数量; d为颗粒粒径; φ为系统中颗粒体积分数, 并假设系统任意颗粒i初始磁矩为m_i.
2 仿真结果与讨论2.1 范德华能对颗粒成链的影响分析随磁性颗粒量级变化, 成链过程中范德华势能对系统相邻颗粒团聚成链的影响可通过微米和纳米级MRFs的Monte Carlo仿真结果进行比较, 对颗粒数64、平均粒径分别为9 nm和2 μm的系统, 在B=0.25 T, φ=10%, 室温条件下, 绘制末状态位形与历次循环能量组成曲线(见图 1, 图 2).
图 1(Fig. 1)

图 1 纳米磁性颗粒系统能量组成与位形Fig.1 Energy composition & nanoparticles' position

图 2(Fig. 2)

图 2 微米磁性颗粒系统能量组成与位形Fig.2 Energy composition & micron particles' position

比较图 1和图 2右侧位形图, 施加外场后的微米和纳米级系统均为多链、单链、团聚等多种情况共存, 颗粒磁矩与外场方向随循环次数增加逐渐趋同, 这与光学显微镜观测结果基本一致^[14].
观察图 1和图 2左侧能量组成曲线发现:不同粒径量级系统按式(4)计算末状态颗粒范德华能, 在总势能中占比差别很大:d=9 nm时, E^vdw=6.64×10^-18 J与极化能E^dd=8.06×10^-18 J在颗粒总能量中的占比接近, 按式(1)计算热耦合系数为15.202 7;d=2 μm时, 比例下降到E_ij^vdw/E_ij^dd1.3×10^-8, 图 2中范德华能为近零直线, 热耦合系数为1.559 9×10⁸, 即:出现极强的成链团聚趋势, 热运动对团聚影响已可忽略不计.
2.2 热耦合系数的试验设计步骤式(1)表达了单独颗粒对之间的磁性能与热能之比, 但因为磁性流体系统结构的复杂性, 系统级热耦合系数为分布粒径d、外磁场强度B以及体积分数φ等参数综合影响下的结果, 单独颗粒对之间的λ难以反映系统级别的能量组成情况, 同时试验也难以观测和采集足够的λ_ij样本集.
故考虑采用如下步骤解决上述问题:
1) 颗粒位形受外磁场强度、颗粒粒径和体积分数三个因子影响, 采用中心组合设计法, 以系统热耦合系数 <λ>为设计响应值, 选择粒径d([6~12 nm]), 体积分数φ([4%~16%])及外场强度B([0.11~0.34 T])作为独立参量, 在不同试验水平组合变量完成DOE方案, 并按照响应面法分析3个主项因子和交互因子对响应值的影响贡献度;
2) Monte Carlo仿真计算不同水平设计参数位形坐标和系统热耦合系数 <λ>值, <λ>由式(8)获得.

(8)

式中的〈E_ij^dd〉是每个颗粒与其他颗粒极化能之和除以单元系统颗粒总数的极化能均值.式(8)与式(1)的区别在于能够考虑单元系统内颗粒间距、局部有道磁场强度等所有因素对极化能的影响.
CCD设计生成的3因子5水平编码及实际值, 如表 1所示.
表 1(Table 1)

表 1 5水平编码及实际响应值Table 1 Five coded levels & actual response values 编号 水平 d/nm φ/% B/T 1 -1.5 6 4 0.11 2 -1 7 6 0.15 3 0 9 10 0.22 4 1 11 14 0.30 5 1.5 12 16 0.34

表 1 5水平编码及实际响应值 Table 1 Five coded levels & actual response values

CCD设计在星点和中心点处共完成2³+6+2=16次试验.用Design Expert^?软件做3因子5水平CCD设计，如表 2所示, 因子已编码且输入仿真试验数据.
表 2(Table 2)

表 2 热耦合系数试验设计与响应值Table 2 Design of experiment & response value of thermal coupling factors 试验轮次编号 试验设计编号 实际值 响应值<λ > d/nm φ/% B/T 1 7 7 14 0.30 0.693 5 2 2 9 16 0.22 2.228 9 3 15 9 10 0.22 2.021 8 4 9 6 10 0.22 0.206 4 5 14 9 10 0.34 1.954 6 6 1 7 06 0.15 0.356 0 7 4 11 14 0.15 2.586 9 8 10 12 10 0.22 2.317 8 9 3 7 14 0.15 0.411 3 10 13 9 10 0.11 1.092 6 11 11 9 04 0.22 0.847 4 12 16 9 10 0.22 1.953 9 13 8 11 14 0.30 3.485 3 14 6 11 06 0.30 2.090 9 15 5 7 06 0.30 0.324 6 16 2 11 06 0.15 1.104 4

表 2 热耦合系数试验设计与响应值 Table 2 Design of experiment & response value of thermal coupling factors

3) 通过方差分析, 缩减无显著回归特征的交互项, 经缩减模型分析的结果如表 3所示.
表 3(Table 3)

表 3 热耦合系数方差检验结果Table 3 ANOVA of thermal coupling factors 因子项 F检验数值 假设概率Prob > F 贡献度 % Model 139.30 <0.000 1? A-d 196.38 <0.000 1? 28.83 B-φ 203.50 <0.000 1? 29.88 C-B 82.86 <0.000 1? 12.16 AB 66.25 0.000 2? 9.73 AC 29.41 0.001 6? 4.32 A2 50.53 0.000 4? 7.42 B2 15.79 0.007 3? 2.32 C2 17.13 0.006 1? 2.51 AB2 13.59 0.010 2? 2.00 失拟度 5.71 0.307 1?? 0.83

表 3 热耦合系数方差检验结果 Table 3 ANOVA of thermal coupling factors

表 3第1列为所有因子主项和交互项，其中A，B和C在方差分析中依次表示直径、体积分数和外磁场强度; 第2列F检验值用于被解释变量在回归方程中的显著性推断; 第3列为显著性水平; 最后一列贡献度是每个因子的F值在全部因子F值总和中的占比, 贡献度越大, 表明该因子在回归模型中对 <λ>值影响越显著.带“?”上标的主因子和交互因子对回归模型影响较显著(高于F值概率大于显著性水平5%); 带“??”上标的最后一行的模型失拟度, 结果30.71%(>10%)说明失拟不显著; 同时模型假设概率是Prob值高于显著水平(5%)的概率, 因此模型对所列因子存在显著的回归特性, 对因子实际值拟合模型见式(9).

(9)

式(9)中的拟合系数A和对应因子项如下:
图 3是按CCD试验设计次序获得的仿真数据和用式(9)计算的16个拟合数据比较, 可以看出二者差异很小.
图 3(Fig. 3)

图 3 热耦合系数拟合与仿真数据比较Fig.3 Comparison of thermal coupling factor's fitting and simulation data

d, φ, B三个变量对热耦合系数λ都有明显影响, 粒径d和体积分数φ分别有28.83%和29.88%的贡献度, 外磁场强度B则只有12.16%, 这证实了之前的设想, 即:颗粒极化能E_ij^dd正比于粒径的立方, 纳米颗粒过小的粒径致使施加外磁场后, 相互间的极化能在数值上与颗粒热运动能相比并不占优势, 当粒径尺度达到≈10 nm时, 增加外磁场B对团聚的促进作用已不再非常明显.此外, 颗粒体积分数具有最高的热耦合系数贡献度, 则存在两个方面的原因:首先, 颗粒相对比例的增加导致颗粒之间局部诱导磁场相互作用区域的叠加, 系统极化能相应增加；此外, 体积分数增加对颗粒自由热运动也是一种几何条件上的限制, 而这一点无法通过式(1)得到体现.
2.3 仿真结果验证Kruse根据一定试验数据, 提出估算成链极限直径的经验方法(式(10))^[8], 计算出室温条件磁流变液在外磁场下达到磁化饱和强度(M_d=0.402 5 T)所生成链状结构的最小颗粒极限粒径约5.7 nm.

(10)

而拟合表达式(9), 限定d, φ和B 3个主因子的范围, 寻找令λ=1的特定解空间.利用Design-Expert^?软件中的优化功能, 设定粒径为最小值, 体积分数与外磁场设为CCD设计范围内, 得如表 4所示的3组优化解.
表 4(Table 4)

表 4 极限粒径的优化解Table 4 Optimal solution of critical CIPs' diameter values No. d/nm φ B/T <λ > 1 7.60(min) 0.073(in range) 0.30(in range) 0.995 5 2 7.09(min) 0.12(max) 0.25(in range) 1.002 1 3 7.063(min) 0.11(in range) 0.27(max) 0.999 6

表 4 极限粒径的优化解 Table 4 Optimal solution of critical CIPs' diameter values

表 4列出三种基本条件:最小(min)、CCD设计范围(in range)和最大(max), 按三种条件约束基本因子, 在满足 <λ>≈1的目标下寻优, 结果显示颗粒在系统中能够成链的极限最小粒径约为7.06 nm, 高于式(10)的估算结果5.7 nm.
拟合方程优化解比文献经验公式高出23%, 说明在实际磁性流体中, 发生流变效应所需的实际粒径条件比仅考虑单独颗粒对的磁热能比要更加严格, 其根本原因在于表达式(9)和式(10)考虑的基本条件不同:式(10)是对单独颗粒对极化能与热运动能之比(λ≈1), 没有纳入颗粒本身受外场环境影响的因素, 无法考虑体积分数对热运动的制约作用, 也不能描述局部诱导磁场对颗粒极化能的增幅; CCD拟合模型考虑了粒径和外磁场、颗粒体积分数等主项, 以及交互项因子的影响, 位形信息包含范德华能、空间位阻等约束条件, 因此极限粒径计算方法能全面反映热运动对磁流变系统整体的影响.
3 结论1) 磁流变系统内热运动对颗粒能量组成的影响程度随粒径的减小而呈立方幂律增大, 纳米量级磁流变液需考虑热运动对极化成链的扰动作用;
2) 粒径、体积分数和外磁场强度等参数间存在交互影响项, 因此单独颗粒对之间的热耦合系数不能准确反映出系统整体的磁、热能量构成情况;
3) 通过Monte Carlo仿真结果中获取颗粒位形信息, 并对外磁场强度、颗粒粒径和颗粒体积分数3个因子在不同水平处进行中心组合设计, 发现外场磁感应强度小于0.3 T、颗粒体积分数小于16%时, 如果颗粒平均粒径小于7 nm, 热运动在颗粒能量组成中占有较大比例, 系统在外磁场作用下将难以形成有效磁链, 宏观上表现为磁致剪力过低, 不构成流变效应发生的基本条件.
参考文献

[1]	Phulé P P. Magnetorheological (MR) fluids:principles and applications[J].Smart Materials Bulletin, 2001, 2001(2): 7–10.DOI:10.1016/S1471-3918(01)80040-X
[2]	Iglesias G R, López-lópez M T, Durán J D, et al. Dynamic characterization of extremely bidisperse magnetorheological fluids[J].Journal of Colloid & Interface Science, 2012, 377(1): 153–159.
[3]	李海涛, 彭向和. 磁流变液铁磁颗粒的动力学分析[J].重庆大学学报, 2010, 33(5): 100–104, 113. ( Li Hai-tao, Peng Xiang-he. Dynamic analysis on ferromagnetic particles of magnetorheological fluids[J].Journal of Chongqing University, 2010, 33(5): 100–104, 113.DOI:10.11835/j.issn.1000-582X.2010.05.016)
[4]	Kormann C, Laun H M, Richter H J. MR fluids with nano-sized magnetic particles technology[J].International Journal of Modern Physics B, 1996, 10(23/24): 3167–3172.
[5]	Heine M C, de Vicente J, Klingenberg D J. Thermal transport in sheared electro-and magnetorheological fluids[J].Physics of Fluids, 2006, 18(2): 023301.DOI:10.1063/1.2171442
[6]	Segovia-Gutiérrez J P, Vicente J D, Hidalgo-álvarez R, et al. Brownian dynamics simulations in magnetorheology and comparison with experiments[J].Soft Matter, 2013, 9(29): 6970–6977.DOI:10.1039/c3sm00137g
[7]	Charles S W. Ferrofluids:magnetically controllable fluids and their applications[M]. Berlin: Springer, 2002.
[8]	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
[9]	Liang M, Wanli S, Rensheng W, et al. Study on shear stress model of magnetorheological fluids with distance weighted factors[J].Smart Materials and Structures, 2017, 26(6): 065009.DOI:10.1088/1361-665X/aa6b5e
[10]	Davis S W, Mccausland W, Mcgahagan H C, et al. Cluster-based Monte Carlo simulation of ferrofluids[J].Physical Review E, 1999, 59(2): 2424–2428.DOI:10.1103/PhysRevE.59.2424
[11]	Donselaar L N, Frederik M P, Bomans P, et al. Visualisation of particle association in magnetic fluids in zero-field[J].Journal of Magnetism and Magnetic Materials, 1999, 201(1/2/3): 58–61.
[12]	Metropolis N, Rosenbluth A W, Rosenbluth M N, et al. Equation of state calculations by fast computing machines[J].The Journal of Chemical Physics, 1953, 21(6): 1087–1092.DOI:10.1063/1.1699114
[13]	Rü?i?ka ?, Allen M P. Collective translational and rotational Monte Carlo moves for attractive particles[J].Physical Review E, 2014, 89(3): 033307.DOI:10.1103/PhysRevE.89.033307
[14]	Jeon J, Koo S. Viscosity and dispersion state of magnetic suspensions[J].Journal of Magnetism and Magnetic Materials, 2012, 324(4): 424–429.DOI:10.1016/j.jmmm.2011.08.025