磁流变液剪切特性的分子动力学模拟

磁流变液剪切特性的分子动力学模拟
王昭轩;丁建国
【摘要】磁流变液是一种高效可控的新型材料,研究其磁流变机理有着广阔的应用前景.目前学者对磁流变剪切特性的多因素分析研究较少.采用LAMMPS分子动力学的方法建立了磁流变液微观模型,分析了磁性颗粒间的相互作用,对磁流变液的剪切过程进行了二维状态下的微观数值模拟,并综合分析了影响剪切屈服应力的主要因素.模拟结果表明,磁性颗粒会沿着磁场方向形成若干条长链状结构,颗粒链受到剪切作用发生变形、倾斜甚至断裂;剪切速率较低时,剪切屈服应力保持稳定;随着颗粒体积百分率的增加,剪切屈服应力呈线性增加的趋势;低磁场强度时,剪切屈服应力随磁场强度增加而增长较快,较高磁场强度时,剪切屈服应力随磁场强度增加的增速放缓;在一定范围内颗粒半径越大,粒径分布越集中,剪切屈服应力越大;随着模拟温度的增加,剪切屈服强度先保持稳定,后下降逐渐明显.
【期刊名称】《功能材料》
【年(卷),期】2019(050)006
【总页数】6页(P46-51)
【关键词】磁流变液;分子动力学;剪切屈服应力
【作者】王昭轩;丁建国
【作者单位】南京理工大学理学院,南京 210094;南京理工大学理学院,南京210094
【正文语种】中文
【中图分类】TB303
0 引言
磁流变液是一种性能可控，反应迅速的新型智能材料。

未施加外部磁场作用时，磁流变液中的颗粒呈随机分布；而在磁场作用下，磁性颗粒会沿着磁场方向呈长链状排列，限制了磁流变液的流动。

目前许多学者借助实验手段来研究磁流变液的机理，Jolly等[1]发现了磁流变液在外加磁场作用下呈链状排列，Liu 等[2]观测了二维情况下不同磁场强度对磁流变液微结构的影响，Furst等[3]观测了颗粒链中的缺陷在链变形过程中的影响，赵春伟[4]观测了不同磁场强度对磁流变液链束形态和间距的影响，但是采用实验的方法难以系统分析磁流变液的演化过程和内在力学关系。

部分学者提出了基于微结构的分析模型。

Ginder等[5]基于链化模型，采用有限元法计算了颗粒间的作用力和磁流变液的剪切屈服应力,Bossi[6]考虑了颗粒饱和磁化和链状结构变形对磁流变液剪切屈服应力的影响，金昀等[7]建立了一种截断范围为圆柱体，磁性颗粒为球体的单元体模型，Shulman等[8]假设磁性颗粒链上的磁矩和外力矩平衡，拟合了磁流变液剪切屈服应力公式，Jolly等[9]考虑了磁性颗粒的极化，计算了颗粒间的自由能和磁流变液剪切屈服应力。

由于磁流变液链化机理复杂，目前借助微结构分析模型来计算颗粒间的剪切屈服应力，难以综合考虑众多因素的影响，难以保证计算精度。

采用分子动力学的方法进行微观数值模拟，可以模拟更符合实际的微结构，考虑各种因素的影响。

杨仕清等[10]考虑了磁流变液中的作用势，利用Metroplis算法模拟了磁流变液的成链模型，Stephen[11]采用不同粒径组分来描述颗粒半径大小，研究了磁流变液链状结构与流动特性，李海涛，彭向和等[12]等建立了微观颗粒的动力学模型，模拟了二维状态下磁性颗粒成链与剪切流动的过程，分析了剪切应变速率的影响。

李杰如，杜成斌等[13]采用微观数值模拟的方法，分析了磁流变
液颗粒百分率的影响。

但目前的微观数值模拟，对于剪切模式下的多因素综合分析研究较少。

本文采用分子动力学软件LAMMPS，模拟了微观条件下磁流变液的剪切模型，综合考虑了剪切应变速率、颗粒百分率、磁场强度、颗粒半径、温度等因素对磁流变液剪切屈服应力的影响，为全面了解磁流变工作机理提供了有效的途径。

1 磁性颗粒受力分析
1.1 颗粒之间的磁作用力
将磁流变液所置的磁场区域视为匀强磁场，将磁性固体颗粒简化为圆球状，质量分布均匀的磁偶极子，颗粒半径为R，体积为V。

假设磁流变载液不可被磁化，忽略磁性固体颗粒被磁化后对外加磁场的影响。

那么颗粒被磁化后的磁矩为
(1)
式中，R为颗粒半径，χ为磁性颗粒的磁化率，H为外加磁场强度。

根据磁偶极子理论，磁性固体颗粒i所受到的磁作用力为
(2)
式中，μ0为真空磁导率，rij为两个磁性固体颗粒i和j之间的相对位置矢量，方向指向i，rij=‖rij‖，mir和mjr分别表示mi和mi沿着相对位置矢量rij方向分量的大小[14]。

假设磁场方向沿y方向竖直向上，两磁性固体颗粒i和j之间的中心连线与磁场方向的夹角为θij，如图1所示。

图1 两颗粒间磁作用力Fig 1 Magnetic force between particles
则磁作用力公式可表示为
(3)
式中，为两磁性固体颗粒相对位置rij的单位矢量，为竖直向上y方向的单位矢量。

1.2 颗粒之间的排斥力
磁性固体颗粒在运动过程中由于相互靠近或碰撞发生排斥，此过程比较复杂，此处采用较为简洁的指数形式来表示颗粒间的排斥力
(4)
式中，β为材料参数，F0取两颗粒接触时的最大磁力，则F0取值为
(5)
1.3 基液对颗粒的粘性阻力
忽略颗粒的转动作用，颗粒在基液中运动受到的阻力由斯托克斯阻力公式表示
(6)
式中，η为磁流变液的粘性系数，v为磁性固体颗粒运动的相对速度矢量。

1.4 颗粒运动的动力学表达式
磁性固体颗粒的动力学表达式为
(7)
2 分子动力学建模
2.1 LAMMPS计算模型的建立
本文使用LAMMPS软件进行分子动力学建模，LAMMPS软件可以模拟多种类型
的分子动力学问题，提供多种力函数和势函数，可以模拟各种系综下的微观颗粒体
系，并进行相应的计算，有较高的计算效率[15]。

本文采用sphere/dipole类型
的磁偶极子，来模拟磁性羰基体粉颗粒，颗粒为圆球状，质量分布均匀[16]。

由于模拟对象颗粒数值单位较小，为了减少计算中的截断误差，对于各个物理量均采用lj类型的约化单位[17]采用磁性固体颗粒的距离参数σ，能量参数ε，质量参
数mi作为基本约化单位，其它物理量均以此为基准进行约化
添加*的量表示约化后的物理量，未添加*的量表示约化前的物理量。

模型x方向采用周期性边界条件，y方向采用固定边界条件，z方向采用周期性边界条件。

外加磁场方向沿y方向，磁流变液流动沿x方向。

初始时刻，颗粒按sq2晶格均匀排布，颗粒间距离按照相应密度给出，颗粒初始
速度和方向根据温度按高斯分布随机确定，整个体统处于微正则系综(NVE)中，采用Verlet积分算法，该方法需要记录每个时间步的速度、受力和坐标，有着较好
的精度和稳定性。

磁性颗粒之间的相互作用采用邻域列表算法，截断半径取2.5R。

当颗粒间距大于
所设截断半径时，默认忽略颗粒间相互作用。

初始NVE积分时，剪应变速率设为0，直到系统达到平衡态。

2.2 LAMMPS计算模型的验证
李海涛[18]等建立了二维状态下磁流变液的微观动力学模型，本文选取其模拟参数，取计算仿真区域磁性固体颗粒半径R=30 μm，颗粒密度ρ=7.5×103 kg/m3，颗粒磁化率χ=1，材料参数β=9，基液的粘度ηj=10-3 Pa·s，真空磁导率
μ0=4π×10-7 H/m，外加磁场强度H=10.7 kA/m，步长Δt=0.1 μs，基于LAMMPS分子动力学的建模方法，模拟了两极板间磁性颗粒成链及剪切流动的过程。

当模型达到稳定状态时，磁流体流动方向与外加磁场方向垂直，磁性固体颗粒沿磁
场方向排列成长链状结构，从而增大了极板间的流动阻力。

在磁性固体颗粒成链的过程中，系统的总磁势能不断减小并趋于固定值，系统微观成链结构逐渐趋于稳定[19]，如图2所示。

同时，模型中除了长通链之外，还存在部分的支链和断链，与李海涛模拟结论相符，如图3所示。

图2 Lammps的链化模拟Fig 2 Chain simulation of Lammps
图3 李海涛的链化模拟Fig 3 Chain simulation of Li
假设模型下极板固定，上极板在外力作用下作剪切运动，剪切应变速率为模拟磁流变液的剪切流动过程，根据LAMMPS分子动力学的方法求出剪切屈服应力随剪切应变速率的变化关系，图4为本文模拟结果与李海涛[18]结论的对比，两者吻合较好。

验证了LAMMPS软件数值模拟可以有效地模拟磁流变液的微观结构。

图4 模拟结果对比Fig 4 Comparison of simulation results
3 磁流变液剪切特性分析
模拟过程中，取计算仿真区域磁性固体颗粒半径R=5 μm，颗粒密度ρ=7.5×103 kg/m3，颗粒体积百分率φ=10%，材料参数β=10，基液的粘度ηj=2×10-3 Pa·s，真空磁导率μ0=4π×10-7 H/m，外加磁场强度H=20 kA/m，温度T=298 K，步长Δt=0.5 μs。

假设模型下极板固定，上极板在外力作用下作剪切运动，剪切应变速率为磁流变液发生剪切流动，如图5所示。

在剪切作用下，磁性颗粒链沿剪切方向倾斜，发生伸长变形。

随着变形的增加，颗粒的间隙逐渐增大，颗粒间的作用力也随之变化，最终磁性颗粒链在作用力最小处发生断裂。

图5 Lammps的剪切模拟Fig 5 Shear simulation of Lammps
3.1 剪切应变速率的影响
取颗粒体积百分率φ=10%，颗粒半径R=5 μm，分别取磁场强度B1=0.2 T，
B2=0.4 T，B3=0.6 T。

在较低的剪切应变速率(10～100 s-1)的情况下，模拟磁流变液的剪切运动，剪切应变速率与剪切屈服应力τ的关系如图6所示。

图6 剪切应变速率与剪切屈服应力Fig 6 Shear strain rate and shear yield stress
采用自行研制的磁流变阻尼器进行测试，其结果与之基本吻合[20]。

从图6可以看出，在匀强磁场的作用下，磁流变液剪切屈服应力与剪切应变速率的关系可以用简单的Bingham黏塑性模型来描述。

磁性颗粒通链受到剪切作用，发生变形甚至断裂。

在剪切应变速率较小时，旧的颗粒链断裂后，磁性颗粒可迅速在截断半径内再次组成新链，磁性颗粒再次成链的速度和断链的速度保持基本平衡，此时随着剪切应变速率的增加，磁流变液微观结构保持稳定，剪切屈服应力的大小保持基本不变。

相同剪切应变速率下，磁场强度越大，剪切屈服应力越大。

其趋势符合Carlson、Catanzarete和St Clair[21]的实验结果。

3.2 磁场强度的影响
根据磁偶极子的理论，磁场强度是影响颗粒间磁作用力的重要因素。

取颗粒半径
R=5 μm，剪切应变速率分别取颗粒体积百分率φ1=10%，φ2=20%，φ3=30%。

在不同的磁场强度的情况下，模拟磁流变液的剪切运动，磁场强度H与剪切屈服
应力τ的关系如图7所示。

采用自行研制的磁流变阻尼器进行测试，其结果趋势与之基本吻合[20]。

从图7可以看出，在磁场强度较小的情况下，磁性颗粒尚未达到磁饱和状态，此时颗粒间的磁作用力较小，表现为磁流变液的剪切屈服应力较小，随着磁场强度的增加，颗粒间的磁作用力增长明显，剪切强度增长较快。

在相对较高的磁场强度时，磁性颗粒饱和达到一定程度，磁场强度的变化对剪切屈服应力的影响变小，大致在B=0.7 T 时，剪切屈服应力增长明显放缓。

磁场强度相同的情况下，颗粒体积百分率增加，剪切屈服应力相应增大。

其趋势符合Jason P Rich等[22]与Daniela Susan-Resiga等[23]的实验结果。

图7 磁场强度与剪切屈服应力Fig 7 Magnetic field intensity and shear yield
stress
3.3 颗粒体积百分率的影响
取颗粒半径R=5 μm，剪切应变速率分别取磁场强度B1=0.2 T，B2=0.4 T，
B3=0.6 T。

在不同颗粒体积百分率的情况下，模拟磁流变液的剪切运动，颗粒体
积百分率φ与剪切屈服应力τ的关系如图8所示。

图8 颗粒体积百分率与剪切屈服应力Fig 8 Particle volume fraction and shear yield stress
从图8可以看出，在颗粒体积百分率取5%～30%的区间内，剪切屈服应力τ随颗粒体积百分率φ的增大而呈线性增加趋势。

颗粒体积百分率取5%时，剪切模型中沿着磁场方向的通链数量较少，主要是未贯通的短链，此时磁流变液的剪切屈服应力较小。

当磁性颗粒体积百分率φ增大时，在外加磁场的作用下，沿磁场方向的
磁通链数目随之增多，相邻磁通链的间距变小，断链和支链比例明显减少，抵抗剪切变形的能力也相应增加，从而提高了磁流变液的剪切屈服应力。

在相同的颗粒体积百分率的情况下，磁场强度增加，剪切屈服应力随之增大。

其趋势符合Jason P Rich等[22]与Daniela Susan-Resiga等[23]的实验结果。

3.4 颗粒半径的影响
取剪切应变速率取颗粒体积百分率φ=20%。

分别取颗粒平均半径R1=5 μm，
R2=5 μm，R3=5 μm。

在每种平均粒径下，分别取不同的粒径组分，服从的正态分布。

在不同的磁场强度的情况下，模拟磁流变液的剪切运动，粒径与剪切屈服应力τ的关系如图9所示。

图9 粒径与剪切屈服应力Fig 9 Particle size and shear yield stress
从图9可以看出，在固体颗粒半径较小时(10 μm以下)，颗粒大小对磁流变液剪
切屈服强度有明显的影响。

在平均粒径从2～10 μm的增加过程中，剪切屈服强
度数值增大了若干倍。

由于半径增加，布朗运动减缓，剪切屈服强度增速逐渐放缓。

粒径分散系对于剪切屈服强度也有一定的影响，平均粒径相同时，正态分布σ2值越小，即半径分布越集中，对应的剪切屈服强度越高，随着磁场强度的增加，该效果越明显。

其趋势符合Lemaire等[24]的实验结果。

3.5 温度的影响
取颗粒半径R=5 μm，剪切应变速率颗粒体积百分率φ1=20%，分别取磁场强度B1=0.4 T，B2=0.5 T，B3=0.6 T。

在不同的温度下，模拟磁流变液的剪切运动，温度与剪切屈服应力τ的关系如图10所示。

图10 温度与剪切屈服应力Fig 10 Temperature and shear yield stress
从图10可以看出，随着温度的上升，磁流变液剪切屈服强度整体呈下降趋势。

当温度20～80 ℃区间时，剪切屈服强度数值基本保持不变；当温度80～140 ℃区间时，剪切屈服强度呈缓慢下降趋势；当温度超过140 ℃时，剪切屈服强度下降趋于明显。

在温度上升的过程中，颗粒的分子布朗运动增强，颗粒的成链效果受到削弱。

取模拟区域统计20和150 ℃情况下的该区域的颗粒数量，后者颗粒数量下降约17%，可知由于热膨胀效应导致模拟间隙内相同体积下的颗粒数量减少，从而使得剪切屈服强度降低。

其趋势符合陈飞等[25]的实验结果。

4 结语
采用分子动力学软件LAMMPS，模拟了微观条件下磁流变液的链化与剪切模型。

根据模拟结果得到了以下结论：
(1) 通过与文献[18]参数进行建模比较及分析，验证了本文采用的LAMMPS分子动力学数值模拟方法，可以有效地模拟磁流变液微观结构。

在匀强磁场作用下，磁性颗粒会沿着磁场方向形成若干条颗粒通链，并存在部分支链和断链。

在剪切作用下，颗粒链沿着剪切运动的方向倾斜，发生伸长变形。

(2) 根据LAMMPS数值模拟，剪切应变速率较低时，磁性颗粒再次成链的速度和断链的速度基本保持平衡，磁流变液微结构保持稳定，剪切屈服应力大小基本保持
不变。

(3) 根据LAMMPS数值模拟，在磁场强度较低时，磁性颗粒磁饱和率相应较低，
颗粒间的磁作用力较小，表现为磁流变液的剪切屈服应力较小。

随着磁场强度的增加，颗粒间的磁作用力增长明显，剪切强度增长较快。

磁场强度相对较高时，磁性颗粒饱和程度较高，此时随磁场强度增加，剪切强度增长缓慢。

(4) 根据LAMMPS数值模拟，在颗粒体积百分率取5%～30%的区间内，剪切屈
服应力τ随颗粒体积百分率φ的增大而呈线性增加趋势。

颗粒体积百分率较低时，沿着磁场方向的通链数量较少，当磁性颗粒体积百分率φ增大时，磁通链数目随
之增多，相邻磁通链的间距变小，断链和支链比例明显减少，抗剪切变形的能力增加，提高了剪切屈服应力。

(5) 根据LAMMPS数值模拟，当粒径处于2～10 μm典型值范围内时，随着固体颗粒半径的增加，颗粒间的磁力增大，磁流变液剪切屈服强度增长明显。

由于半径增加，布朗运动效应减小，颗粒沉淀趋势增大，导致了剪切屈服强度增速逐渐放缓。

粒径分散系对于剪切屈服强度也有一定的影响。

平均粒径相同时，粒径正态分布
σ2越小，即粒径分布越集中，剪切屈服强度越高。

(6) 根据LAMMPS数值模拟，在温度低于80 ℃时，磁流变液剪切屈服强度保持
基本稳定，在温度超过80 ℃时，由于磁流变液热膨胀现象，单位体积的磁流变颗粒减少，并且颗粒受热磁化率降低，从而导致了剪切屈服强度下降逐渐明显。

参考文献：
【相关文献】
[1] Jolly M R, Carlson J D, Munoz B C.A model of the magnetorheological
materials[J].Smart Materials & Structures, 1996, 5(5):607-614.
[2] Liu X H, Liu Y Z, Liu H.Magnetorheological fluids particles simulation through integration of Monte Carlo method and Gpu accelerated technology[J].CMES: Computer Modeling in Engineering & Sciences, 2013, 91(1):65-80.
[3] Furst E M, Gast A P.Dynamics and lateral interactions of dipolar chains[J].Physical Review E, 2000, 62(5):6916-6925.
[4] Zhao Chunwei. Study on mechanical properties of magnetorheological fluids based on microstructure[D].Chongqing: Chongqing University, 2014:1-115(in Chinese).
赵春伟.基于微结构的磁流变液力学性能研究[D]. 重庆: 重庆大学, 2014:1-115.
[5] Ginder J M, Davis L C, Elie L D.Rheology of magnetorheological fluids: models and measurements[J].International Journal of Modern Physics B, 1996, 10(23):3293-3303. [6] Bossis G, Lacis S, Meunier A, et al.Magnetorheological fluids[J].Journal of Magnetism and Magnetic Materials, 2002, 252(1):224-228.
[7] Jin Yun, Zhang Peiqiang, Wang Xiaohua, et al.Numerical calculation of shear yield stress of magnetorheological fluid [J].Journal of China University of Science and Technology, 2002, 252(1):224-228(in Chinese).
金昀, 张培强, 汪小华, 等.磁流变液剪切屈服应力的数值计算[J].中国科学技术大学学报, 2001,
31(2):168-173.
[8] Shulman Z, Kordonsky V, Zaltsgendler E, et al.Structure, physical properties and dynamics of magnetorheological suspensions[J].International Journal of Multiphase Flow, 1986, 12(6):935-955.
[9] Jolly M R, Carlson J D, Munoz B C.A model of the magnetorheological
materials[J].Smart Materials & Structures, 1996, 5(5):607-614.
[10] Yang Shiqing, Zhang Wanli, Gong Jie, et al.Study on rheological properties of magnetorheological fluids[J].Journal of Functional Materials, 1998, 10(5):550-552(in Chinese).
杨仕清, 张万里, 龚捷, 等.磁流变液的流变学性质研究[J].功能材料, 1998, 10(5):550-552. [11] Sherman S G, Wereley N M.Effect of particle size distribution on chain structures in magnetorheological fluids[J].Magnetics, IEEE Transactions, 2013, 49(7):3430-3433. [12] Li Haitao, Peng Xianghe, Chen Weimin.Shear yield stress model of magnetorheological fluids based on chain analysis[J].Journal of Chemical Physics, 2005, 18(4):505-509(in Chinese).
李海涛, 彭向和, 陈伟民.基于链化分析的磁流变液剪切屈服应力模型[J].化学物理学报, 2005,
18(4):505-509.
[13] Li Jieru, Du Chengbin. Numerical simulation of microstructures of magnetorheological fluids[J].Journal of Hohai University, 2011, 39(3):311-316(in Chinese). 李杰如, 杜成斌.磁流变液微观结构的数值模拟[J].河海大学学报, 2011, 39(3):311-316.
[14] Liu Y Z, Liu X H, Liu H.Dynamic simulation for three-dimensional micro-structure of
magnetorheological fluids[J].AISS, 2012, 12(4):55-61.
[15] Barski M, Chwal M, Kedziora P.Molecular dynamics in simulation of magneto-rheological fluids behavior[J].Key Engineering Materials, 2013, 542(1):11-27.
[16] Pei L, Pang H M, Ruan X H, et al. Magnetorheology of a magnetic fluid based on
Fe3O4 immobilized SiO2 core-shell nanospheres: experiments and molecular dynamics simulations[J].Rsc Advances, 2017, 7(1):8142-8150.
[17] Li Jie, Peng Yongbo.Stochastic multiscale modeling of dynamic yielding behavior of magnetorheological fluids[J].Computational Physics, 2012, 29(1):95-100(in Chinese).
李杰, 彭勇波.磁流变液动态屈服行为的随机多尺度建模[J].计算物理, 2012, 29(1):95-100. [18] Li Haotao, Peng Xianghe, Chen Weimin.Numerical simulationn and analysis of rheological characteristics of magnetorheological fluid[J].Journal of Functional Materials, 2006, 5(37):710-715(in Chinese).
李海涛, 彭向和, 陈伟民.磁流变液流变特性的数值模拟分析[J].功能材料, 2006, 5(37):710-715. [19] Yi Chengjian, Peng Xianghe, Li Haitao.Morphological stability analysis of magnetorheological microstructures under static magnetic field[J].Journal of Functional Materials, 2008, 12(39):1961-1964(in Chinese).
易成建, 彭向和, 李海涛.静磁场下磁流变微结构形态稳定性分析[J].功能材料, 2008, 12(39):1961-1964.
[20] Zhang Lifeng. Optimization of PID control algorithms based on MR damping structure[D].Nanjing: Nanjing University of Science and Technology, 2016:1-103(in Chinese)
张立峰.基于磁流变阻尼结构的PID控制算法优化研究[D]. 南京: 南京理工大学, 2016:1-103. [21] Carlson J D, Catanzarete D M, St Clair K A. Commercial magneto-rheological fluid devices//Proceedings 5th Int, Conf on ER Fluids, MR Suspensions and Associated Technology[C]. Singapore: World Scientific, 1996.
[22] Rich J R, Doyle P S, Mckinley G H. Magnetorheology in aging, yield stress matrix fluid[J].Rheologica Acta, 2012, 51(1):579-593.
[23] Susan-Resiga D, Vekas L, Mckinley.Yield stress and flow behavior of concentrated ferrofluid-based magnetorheological fluids: the influence of composition[J].Rheologica Acta, 2014, 53(1):645-653.
[24] Lemaire E, Meunier A, Bossis G, et al.Influence of the particle size on the rheology of magnetorheological fluids[J].Journal of Rheology, 1995, 39(5):1011-1020.
[25] Chen Fei, Tian Zuzhi, Wang Jian.Effect of temperature on rheological properties of magnetorheological fluids[J].Journal of Functional Materials, 2014, 45(20):95-98(in Chinese).
陈飞, 田祖织, 王建.温度对磁流变液流变性能的影响研究[J].功能材料, 2014, 45(20):95-98.。