磁流變液微觀結(jié)構(gòu)的數(shù)值模擬

李杰如，杜成斌

(1.河海大學文天學院，安徽馬鞍山 243031;2.河海大學力學與材料學院，江蘇南京 210098)

磁流變液是一種非常具有發(fā)展前途和工程應用價值的新興智能材料，被認為是最具前途的智能材料之一[1].磁流變液是由非膠體的細小顆粒分散溶于絕緣基液中形成的懸浮液，在外磁場作用下能產(chǎn)生明顯的磁流變效應，其力學、熱學、流變學性能等會迅速發(fā)生變化，響應時間很短，為毫秒量級，且這種液態(tài)和固態(tài)之間快速轉(zhuǎn)化是可逆的.在該轉(zhuǎn)化過程中，磁流變液的黏度保持連續(xù)、無級變化，固化強弱可由外加磁場控制，能耗極小，可實現(xiàn)實時的主動控制[2-5].

1948年Rabinow[6]最早發(fā)明了磁流變液及其應用裝置離合器.Ginder等[7]認為在外磁場作用下，磁流變液中的鐵磁性固體顆粒形成的是一條一條的單鏈，并從磁性固體顆粒的非線性磁化、磁飽和以及相鄰顆粒接觸位置處的局部磁場等方面進行了有限元計算，得出了磁流變液的剪切屈服應力的極值與磁性固體顆粒的飽和磁化強度成正比的結(jié)論，該結(jié)論也被后來的實驗所證實.Rosensweig[8]設計出一種平均場模型——層模型，由于該模型忽略了相鄰顆粒間的磁場集中效應，其結(jié)果與實驗值相比往往偏低.Bossis等[9-10]將Klingenberg等[11]利用電流變液推導出的恢復力公式應用于磁流變液，但由于沒有考慮磁性粒子的非線性磁化過程，其結(jié)果往往偏大.

近年來，國內(nèi)武漢理工大學的阮中尉[12]建立了BCT結(jié)構(gòu)(體心四方結(jié)構(gòu))模型和物理材料模型，得到了磁流變液受小角度剪切變形情況下其剪切應力的理論計算公式;中國礦業(yè)大學的趙四海等[13]研究了磁流變液中顆粒的相互作用，分析了磁流變液產(chǎn)生鏈狀結(jié)構(gòu)的原因，推導出磁流變液在外磁場作用下的屈服應力公式;廣東工業(yè)大學陳偉俊[14]在對各種磁流變液的力學模型對比分析的基礎上，根據(jù)磁場的庫侖定律，提出了一種能夠簡潔表示各種因素對磁流變液力學性能影響的力學模型.

基于微結(jié)構(gòu)的建模及力學分析只能描述磁流變液的宏觀特性，不能描述微結(jié)構(gòu)的形成機制及其演化過程.影響磁流變液性能的因素很多，單純從實驗或宏觀角度很難系統(tǒng)分析各因素的影響及效果，而基于微觀角度的數(shù)值模擬方法可以有效彌補這些缺陷，有利于優(yōu)化磁流變材料的性能，設計高性能磁流變液.本文從數(shù)值模擬角度出發(fā)，利用分子動力學中的速度Verlet算法，基于磁流變液中磁性顆粒的實際受力，建立磁流變液鏈化和剪切的分析模型，采用Matlab語言編程，分析粒子的運動規(guī)律，探討磁流變液微觀結(jié)構(gòu)的鏈化形成過程和剪切變形過程.

1 磁性顆粒的力學模型

在外加磁場的作用下，顆粒被磁化成偶極子在基液中產(chǎn)生運動，粒子在運動過程中起重要作用的力是顆粒之間的磁力、排斥力和基液對顆粒的黏性阻力.其他力由于影響很小，可以忽略不計.

1.1 磁性顆粒受力分析

1.1.1 顆粒之間的磁力

假設磁流變液中的基液是不可磁化的，其中的鐵磁性固體顆粒為大小相同的圓球狀，粒子的半徑均為R，體積為V，在外加磁場H作用下，把固體顆?？闯墒且粋€個磁偶極子，其方向與外加磁場方向相同.則顆粒被磁化后的磁矩[15]為

式中:M——磁化強度，M=χ H;χ——顆粒的磁化率.

顆粒i受到的磁力表達式[16]為

式中:μ0——真空磁導率;rij——兩顆粒i和j之間的相對位置矢量(指向i)，rij=‖rij‖;mir，mjr——磁矩mi和mj沿相對位置矢量rij方向的分量大小.

假設磁場方向為豎直向上，兩顆粒i和j之間的中心連線與磁場的夾角為θij，利用式(1)，式(2)所示的磁力表達式可改寫為

1.1.2 顆粒之間的排斥力

顆粒在運動過程中由于發(fā)生碰撞而相互排斥，排斥力[17]為

式中β為材料參數(shù).

由式(3)可以看出，當粒子沿磁場方向排列且接觸時，顆粒之間的磁力最大.F0可以取為兩個顆粒相接觸時的最大磁力，此時粒子沿磁場方向排列且接觸，rij=2R，F(xiàn)m=-Fr，F(xiàn)0的取值為

1.1.3 基液對顆粒的黏性阻力

假設顆粒為半徑為R的圓球形，基液為黏性不可壓流體，速度不太大時顆粒在靜止液體中平動時受到的阻力Fv由Stokes阻力公式[18]來表示:

式中:η——液體的黏性系數(shù);r——顆粒的位置矢量;v——顆粒運動的速度矢量.

1.2 顆粒運動的動力學方程

設單個顆粒質(zhì)量為m，則

式中ρ為鐵磁性固體顆粒材料的密度.

由牛頓第二定律，顆粒運動的動力學方程為

1.3 顆粒運動的模擬迭代算法

無外加磁場時顆粒在基液中隨機分布，受到磁場作用后，顆粒即被磁化成一個個的磁偶極子，受到其他顆粒對它的磁力和排斥力等力的作用，從而具有了加速度，這樣速度的增量和位移的增量就可以求出.而顆粒的位置在位移增量的影響下就發(fā)生了改變，從而顆粒受力也發(fā)生了改變，引起加速度和速度的改變.依次求出每個時間步長內(nèi)的加速度、速度和位移，直到系統(tǒng)達到穩(wěn)定即平衡狀態(tài)為止.

分子動力學中，常見的求解運動方程的算法有Verlet算法、Leapfrog(蛙跳)算法、Gera(預測-校正)算法和Velocity Verlet算法(速度Verlet算法)等.速度Verlet算法可以同時給出位置、速度和加速度，并且不損失精度，所以本文的模擬過程選擇速度Verlet算法.該算法需儲存每個時間步顆粒的坐標、速度和力，每步計算涉及2個時間步，需要計算顆粒的位置更新以后和速度更新之前的力[19].

由于無外加磁場時，顆粒是隨機分布在基液中且?guī)缀蹯o止不動，所以顆粒的初始位置由隨機函數(shù)產(chǎn)生，初始速度為零，則t時刻的顆粒運動的加速度為

設模擬時所用的時間步長為Δt，則t+Δt時刻的顆粒位置為

t+Δt時刻的顆粒速度為

顆粒的位置、速度和受力不斷更新，式(9)，(10)和(11)不斷循環(huán)，直至顆粒的位置達到穩(wěn)定狀態(tài)(穩(wěn)定判據(jù)為系統(tǒng)動能小于設定值)，即認為循環(huán)終止.模擬區(qū)域采用周期性邊界條件.

2 模擬結(jié)果分析

2.1 鏈化過程的模擬

模擬中顆粒磁化率 χ=1，真空磁導率 μ0=4π×10-7H/m，材料參數(shù) β=9，固體顆粒材料的密度ρ=7.5t/m3，基液的黏性系數(shù)η=0.001Pa?s，外加磁場強度H=10.7kA/m.為了得到最好的鏈化效果，取顆粒體積百分率ξ=30%，顆粒半徑R=5μ m，步長20μ s.鏈化過程模擬結(jié)果見圖1.顆粒鏈化過程先快后慢，整個鏈化過程所耗費的時間在2～3ms之間.

圖1 顆粒鏈化過程Fig.1 Illustration of chain process for particles

顆粒的體積百分率對磁流變液的性能有很大的影響，一般要在15%～30%范圍內(nèi).下面以R=5 μ m為例，適當擴大模擬范圍來比較不同顆粒體積百分率時顆粒成鏈的區(qū)別.圖2為顆粒體積百分率ξ=15%，20%，25%和30%時的模擬結(jié)果.

圖2 不同顆粒體積百分率的模擬結(jié)果Fig.2 Simulated results of different volume percentages for particles

從圖2可以看出，顆粒體積百分率較小時，所形成的通鏈較少，甚至幾乎沒有，達到穩(wěn)定狀態(tài)時，大多數(shù)粒子還是以支鏈、孤鏈或單個粒子的形態(tài)存在.隨著顆粒體積百分率的增加，通鏈逐漸增多，支鏈、孤鏈和單個的粒子在減少.顆粒體積百分率越大，所成通鏈越多.所以，在條件容許的范圍內(nèi)，在充分考慮粒子沉降性的基礎上，應盡可能地增大顆粒的體積百分率.

2.2 剪切過程的模擬

磁流變液在受到外界的作用發(fā)生剪切流動時，基液也會同時流動.設其沿x正方向流動，剪切速率為﹒γ，則顆粒所受黏性阻力公式(6)應改寫為

以圖2(d)的模擬結(jié)果作為初始條件，上板在外力作用下作勻速剪切運動，剪切速率為﹒γ=200s-1，下板固定不動，其微觀結(jié)構(gòu)的演變過程如圖3所示.

圖3 低剪切速率下的剪切過程模擬Fig.3 Shear process at low shear velocity

由圖3可以看出鏈逐漸向剪切方向傾斜，部分鏈發(fā)生斷裂，同時生成新鏈，新鏈由斷鏈、支鏈和孤鏈等重組生成，但是支鏈、孤鏈和單個粒子幾乎始終存在.舊鏈的斷裂和新鏈的生成基本上達到動態(tài)平衡的狀態(tài)，也就是說，宏觀上，磁流變液的剪切應力基本上應該是穩(wěn)定的.

仍然以圖2(d)的模擬結(jié)果作為初始條件，剪切速率為﹒γ=600s-1，其微觀結(jié)構(gòu)的演變過程如圖4所示.

圖4 高剪切速率下的剪切過程模擬Fig.4 Shear process at high shear velocity

圖4結(jié)果表明，高剪切速率下，粒子向剪切方向的傾斜程度變大，但是通鏈減少，新鏈的生成速度明顯跟不上舊鏈的斷裂速度，在宏觀上應表現(xiàn)為高剪切速率下，隨著剪切速率的增加，磁流變液的剪切應力下降.

2.3 剪切應力的計算

根據(jù)磁力的計算公式，當顆粒沿磁場方向排列時，底板所受磁力在垂直于磁場的方向，即 x方向上分量為零.當顆粒受到外界的剪切作用而傾斜時，磁力在x方向上就有了分量，顆粒對底板在水平方向上的力的分量大小，宏觀表現(xiàn)為剪切應力的大小.

設不受剪切時同一鏈中相距k個顆粒的2個顆粒之間的距離為 dk=2kR，則剪切過程中，鏈被拉長后，該距離變?yōu)閐k=2kR/cosθ，則單鏈對宏觀剪切應力的貢獻，即單鏈的剪切強度為

單位面積內(nèi)的鏈條數(shù)為

式中:S——磁流變液的底面積;L——磁流變液的高度.

單位橫截面面積內(nèi)各鏈對宏觀剪切應力貢獻的和為

因為磁流變液中的磁鏈并不都是通鏈，所以要乘以一個折減系數(shù) φ.

磁流變液的宏觀剪切應力為無鐵磁性顆粒時基液的剪切應力τ0加上鐵磁性顆粒對剪切應力的貢獻ˉτ，即

利用式(13)計算的單鏈的應力-應變曲線見圖5;根據(jù)式(16)計算出的磁流變液的應力-應變曲線見圖6;不同磁感應強度的應力-應變曲線見圖7;不同顆粒體積百分率剪切應力變化情況見圖8.

圖5 單鏈的應力-應變曲線Fig.5 Stress-strain curve of single chain

圖6 磁流變液的應力-應變曲線Fig.6 Stress-strain curve of MFR

圖7 不同磁感應強度的應力-應變曲線Fig.7 Stress-strain curves under different magnetic induction intensities

圖8 不同顆粒體積百分率的剪切應力變化Fig.8 Variation of shear stress under different volume percentages of particles

從圖5、圖6可以看出，不同于單鏈的剪切強度先增大后減小，磁流變液的剪切強度只增大不減小.從圖7可以看出，隨著磁感應強度的增大，磁流變液的剪切應力有著顯著提高.從圖8可以看出，隨著顆粒體積百分率的增大，磁流變液的剪切應力幾乎呈線性增加.

3 結(jié) 論

a.顆粒鏈化的變化過程先快后慢，在1ms之內(nèi)就可以形成鏈化，然后慢慢調(diào)整位置，最終趨于穩(wěn)定，整個鏈化過程所耗費的時間在2～3ms之間.

b.不同顆粒體積百分率的磁流變液成鏈情況不同，顆粒體積百分率越大，模擬時所成通鏈越多，支鏈、孤鏈和單個粒子越少.因此，在充分考慮粒子沉降性的基礎上，在條件容許的范圍內(nèi)，應該盡可能地增大顆粒的體積百分率.

c.磁流變液受到剪切時，鏈逐漸向剪切方向傾斜，部分鏈發(fā)生斷裂，同時斷鏈、支鏈、孤鏈和單個粒子等又重組生成新鏈，低剪切速率下，舊鏈的斷裂和新鏈的生成基本上達到動態(tài)平衡的狀態(tài)，但是支鏈、孤鏈和單個粒子幾乎始終存在.

d.高剪切速率下，粒子向剪切方向的傾斜程度變大，但是通鏈減少，新鏈的生成速度明顯跟不上舊鏈的斷裂速度，與宏觀表現(xiàn)高剪切速率下，隨著剪切速率的增加，磁流變液的剪切應力下降相吻合.

e.單鏈的剪切強度先增大后減小.磁流變液的剪切屈服應力隨磁感應強度的增大而增大，隨顆粒體積百分率的增大而增大.

[1]汪建曉，孟光.磁流變液研究進展[J].航空學報，2002，23(1):6-12.(WANG Jian-xiao，MENG Guang.Research advances in magnetorheological fluids[J].Acta Aeronautica Et Astronautica Sinica，2002，23(1):6-12.(in Chinese))

[2]汪建曉，孟光.磁流變液裝置及其在機械工程中的應用[J].機械強度，2001，25(1):50-56.(WANG Jian-xiao，MENG Guang. Magnetorheological fluid devices and their applications in mechanical engineering[J].Journal of Mechanical Strength，2001，25(1):50-56.(in Chinese))

[3]李金，張華良.磁流變液研究與應用[J].上海大學學報:自然科學版，2004，10(1):21-25.(LI Jin，ZHANG Hua-liang.Study and application of magnetorheological fluids[J].Journal of Shanghai University:Natural Science Editions，2004，10(1):21-25.(in Chinese))

[4]張平.磁流變液的研發(fā)與應用[J].功能材料信息，2006，3(1):19-22.(ZHANG Ping.Study and application of magnetorheological fluids[J].Functional Materials Information，2006，3(1):19-22.(in Chinese))

[5]裴暢貴，郭張霞.磁流變液特性分析[J].機械管理開發(fā)，2007(2):24-25.(PEI Chang-gui，GUO Zhang-xia.Behavior analysis of magneto-rheological fluids[J].Mechanical Management and Development，2007(2):24-25.(in Chinese))

[6]RABINOW J.The magnrtic fluid clutch[J].ALEE Transaction，1948，67:1308-1315.

[7]GINDER JM，DAVID L D.Shear stress in magnetorheological fluids[J].Role of Magnetic Saturation Appl Phys，1994，65(26):2410-2412.

[8]ROSENSWEIG R E.On magnetorheology and electrorheology as states of unsymmetric stress[J].J Rheol，1995，39(l):179-192.

[9]BOSSIS G，LNAIRE E.Yield stresses in magnetic suspensions[J].J Rheol，1991，35(7):1345-1354.

[10]LEMAIRE E，BOSSIS G.Yield stress andwall effects in magnetic colloidal suspensions[J].J Phys D:ApplPhys，1991(24):1473-1477.

[11]KLINGENBER G D J，ZUKOSKI C F.Studies on the steady-shear behavior of electrorheological suspensions[J].Langmuir，1990(6): 15-24.

[12]阮中尉.基于BCT模型的磁流變液剪切應力預測理論研究[D].武漢:武漢理工大學，2006.

[13]趙四海，方佳雨.一種磁流變流體鏈化模型的研究[J].潤滑與密封，2006(11):108-110.(ZHAO Si-hai，F(xiàn)ANG Jia-yu.Theoretics study on the model of magnetic flocculation of magneto-rheological fluid[J].Lubrication Engineering，2006(11):108-110.(in Chinese))

[14]陳偉俊.適用于直線電機阻尼控制的磁流變液的技術研究[D].廣州:廣東工業(yè)大學，2007.

[15]FUST E，GAST A.Micrlmechanics of magnetorheological suspensions[J].Physical Review E，2000，61(6):6732-6739.

[16]方生，張培強.旋轉(zhuǎn)磁場作用下磁流變液顆粒運動及結(jié)構(gòu)演化的模擬[J].化學物理學報，2001，14(5):562-566.(FANG Sheng，ZHANG Pei-qiang.Simulationof the structure and the dynamicsof the particles ofmr fluids in rotating magnetic fields[J].Chinese Journal of Chemical Physics，2001，14(5):562-566.(in Chinese))

[17]SONIA M，CALDERON S.Rotational dynamics in dipolar colloidal suspendions:video microscopy experiments and simulations results[J]. Non-Newtonian FluidMech，2002，102:135-148.

[18]劉望一.流體力學:下[M].北京:北京大學出版社，2000:46-90.

[19]張躍，谷景華.計算材料學基礎[M].北京:北京航空航天大學出版社，2007:50-89.