修正磁化模型的多組分鐵磁性顆粒運動研究1)

陳巨輝 安 然 舒崚峰 李 丹 劉曉剛 毛 穎 陳紀元 高浩銘 呂文生 孟凡奇

* (哈爾濱理工大學機械動力工程學院,哈爾濱 150080)

? (中國電建集團華東勘測設計研究院有限公司,杭州 311122)

** (中國航發哈爾濱東安發動機有限公司,哈爾濱 150066)

引言

鐵磁性顆粒(簡稱“磁性顆?！?在磁場中會自發地進行磁化,其被廣泛地運用于靶向給藥[1-2]、磁流變液[3-5]、磁性分選[6-9]、磁性新材料[10-11]和磁流化床[12-14]等領域中.常規的多相流動過程中產生的混合不均、溝流、節涌和分層等不良現象[15],可以通過添加外部磁場的方法對顆粒運動進行干預,對磁性顆粒流化過程中產生的不良現象加以改善與控制,減少顆粒運動中的不良現象,提高顆粒循環的運行效率,降低顆?？刂频倪\行成本[16].

磁性顆粒在磁場中受到磁場力與磁化力的共同作用,其中磁化力為國內外學者的研究重點.Rosensweig[17]對于單個球形磁性顆粒在磁場中受到的磁場力進行了推導與實驗,結果表明,球形磁性顆粒在非均勻磁場中受到力的作用,該作用的方向與磁場梯度方向一致,其大小與磁場強度梯度、顆粒磁化率和顆粒體積成正比,進一步定量得到Rosensweig 磁場力公式,該計算磁場力模型具有較高的精度.

鐵磁性顆粒在施加外磁場后可認為瞬間完成磁化過程,由于顆粒在磁化后具有極性,因此成對的鐵磁性顆粒在磁場中會相互吸引或者排斥,最終形成聚鏈.Pinto-Espinoza[18]對于鐵磁性顆粒在磁場中的運動進行了理論推導和模擬研究,其假定完全相同的兩個顆粒,在垂直方向均勻外磁場的條件下,兩個理想偶極鐵磁性顆粒之間,施加的外磁場與鐵磁性顆粒磁場疊加,結合磁場勢能公式,推導提出了P-E磁化公式.為整個鐵磁性顆粒模擬研究奠定了基礎,但由于其采用均勻垂直磁場并未考慮磁感應強度梯度和磁場方向對于顆粒的作用,因此該模型具有一定的局限性.Hao 等[19]運用P-E 磁化模型與離散模型相結合,研究鼓泡流化床中添加垂直方向二維梯度磁場的運動狀態,分析了不同磁感應強度下的鐵磁性顆粒的速度與濃度變化.Ke 等[20]運用P-E 磁化模型與離散模型相結合,研究添加二維外部垂直恒定磁場計算域,模擬流體中磁性粒子的運動,得出結論,磁場的方向和磁性粒子的相對位置決定了相互作用力是吸引還是排斥,隨著距離的增加,磁力迅速減小;當施加外部垂直恒定磁場時,磁場中的顆粒會發生鏈式聚集,鐵磁性顆粒形成鏈狀結構.

磁場流化床的眾多優勢使其在工業生產領域包括物料分離、顆粒篩選及改善黏性顆粒流化效果等方面承擔了重要的作用[21-22].對于施加外磁場的顆粒多相流運動,主要研究磁性顆粒在磁流化床內的運動狀態.Valverde 等[23]通過采取將氣流與施加的外部磁場進行模擬耦合,研究磁場輔助下的顆粒流化狀態,結果表明,施加的磁場使得顆粒與顆粒間產生力的作用,沿著磁場方向形成鏈狀結構,顆粒成鏈長度與磁場強度、粒子磁化率、顆粒大小和顆粒密度有關,并提出了考慮形成磁鏈長度的Richardson-Zaki 方程.Zhu 等[24]針對具有可磁化和不可磁化顆粒的流化床,建立了一種基于快速測量相應沉降床內磁導率的新方法.避免了從沉降床中取樣顆粒和隨后分離兩種類型的顆粒,從而消除了這兩個步驟中所關注的許多問題,并且更加高效.楊慧等[25]研究摻雜鐵粉的SiO2顆粒在磁流化床的運動,探究了磁場強度對床內氣含率和局部壓力波動均方根的影響,結果表明,隨著床內磁場強度的增加,局部氣含率和平均氣含率均有所增大;當磁場強度一定時,顆粒粒徑越大,局部氣含率和平均氣含率均增大;當磁場強度和床高增加時,床內壓力波動升高.劉金平[26]研究納米級鐵磁性顆粒在微小磁流化床內的運動,模擬了納米鐵磁性顆粒在不同的床徑、氣速、配風和磁場強度下的運動,結果表明,在內徑20 mm的微小流化床中,施加的外磁場可以抑制納米顆粒產生的溝流現象,幫助顆粒維持穩定流動狀態.Chen等[27]通過計算顆粒流體力學,研究了鼓泡床密相區顆粒的混合擴散對爐內反應速率的影響.Song 等[28]采用雙流體模型結合多組分傳質模型和氣泡中尺度阻力模型來揭示活性顆粒的傳質特性.Ganzha 等[29]對不同粒徑、料床高度的鐵磁性顆粒進行鼓泡流化實驗,結果表明,當施加磁場的磁感應強度足夠大時,觀測到床內產生顯著的磁場,該磁場影響施加的外部磁場,導致外磁場發生變形,鐵磁性顆粒之間產生明顯的作用力,并出現了成鏈團聚現象.Jovanovic等[30]在施加垂直方向梯度磁場的條件下,對于鐵磁性顆粒的流化過程進行模擬,結果表明,添加垂直含梯度磁場,由于產生的磁場力與顆粒所受重力相互抵消,可以實現在微重力環境下的顆粒運動,并且磁力在阻止鐵磁性顆粒在微重力條件下從床中逸出起著主要作用.

綜上所述,學者們對鐵磁粒子運動的研究多停留在垂直磁場和水平磁場上.本文提出一種改進的P-E 磁化模型,將適用范圍從垂直磁場擴大到任意方向,能夠計算空間中任意磁場方向下顆粒所受到磁化力,使得研究任意磁場方向對于磁性顆粒運動的影響成為可能.基于改進的P-E 磁化模型,采用有限體積法-離散單元法(finite volume method-discrete element method,FVM-DEM)模擬,研究了多組分顆粒在磁場中的運動,分析了不同磁場方向和不同配比對多組分顆粒聚鏈的影響.仿真結果為多組分鐵磁粒子在磁場中的運動提供了理論依據和數據支持.為進一步研究鐵磁性顆粒在不同磁場環境中的運動提供了一定的思路.

1 數學模型

1.1 FVM-DEM 模型

氣相采用有限體積法,通過質量守恒方程和動量守恒方程描述.質量守恒方程

動量守恒方程

式中,εg為氣相體積分數,Kgp為單位體積網格內固相與氣相的相間交換系數.

固相采用離散單元法,通過牛頓定律得到顆粒的運動方程與轉動方程.顆粒運動方程

顆粒轉動方程

其中產生的曳力采用Spherical 曳力模型進行描述,其表達式如下

式中,C為曳力系數,a1,a2和a3為根據不同Re的定值常數.其中液固相間耦合需計算由于曳力作用產生的動量交換源項,其表達式如下

顆粒運動發生碰撞受到接觸力Fc作用,采用Hertz-Mindlin[31-32]接觸模型計算顆粒接觸力

式中,Fr為徑向接觸力,Frd為徑向阻尼力,Ft為切向接觸力,Ftd為切向阻尼力,E*為當量楊氏模量,R*為當量彈性模量,δr為徑向重疊量,δt為法向重疊量,m*為當量質量,vn、vt為相對速度的法向、切向分量,α0為修正恢復系數,Dn和Dt為法向和切向剛度.將作用力與阻尼力矢量疊加進而計算出顆粒碰撞合力如下

式中,r為徑向方向向量,t為切向方向向量.

顆粒所受磁場力Fme與磁場方向保持一致,磁化力Fmi可分解為徑向磁化力Fmi,r和切向磁化力Fmi,φ.因此顆粒所受總磁力表達式如下

式中,Vp為顆粒體積,μ0為真空磁導率,χp為磁化率,H為磁感應強度.

1.2 修正P-E 模型

本文采用相對坐標系轉換方法,在相對坐標系中采用傳統P-E 模型,從而實現P-E 模型適用于任意方向磁場.基于常規坐標系xoy,以計算顆粒為圓心,以磁場方向為y'正方向,x'正方向為y'順時針旋轉90°方向,建立相對參考系x'oy',顆粒對相對坐標系如圖1 所示.

圖1 顆粒對相對坐標系Fig.1 Relative coordinate system of particle pairs

圖1 中,β為磁場偏角,γ為偶極角,B0為磁通量密度.由于施加外加磁場方向β為任意方向,可以通過方向向量方式表示磁化偏角和磁偶極矩關系

同理可以通過單位向量表示顆粒中心連線方向

通過相對參考系計算顆粒對產生的磁偶極矩為

其中空間一點磁偶對產生的磁場大小計算公式如下

結合式(16)、式(17)和式(20),通過三角函數向量計算并整理得到

假設兩種顆粒種類大小完全相同,因此m=m1=m2,并代入式(18)得到磁偶極矩表達式

由于選取相同顆粒,通過式(16)與式(22)方向向量系數對照,得到二元一次方程組

通過求解方程組,可解出相對參考系中γ偶極角與m磁偶極矩

因此就可以通過勢能對徑向和角向進行微分,計算出徑向磁力和角磁力

式中U為偶極子之間的勢能.

進而計算徑向磁化力與切向磁化力大小

式(29)和式(30)即為修正P-E 磁化模型.

2 模擬條件

2.1 仿真模型

模擬結構采用偽二維的幾何結構,該結構長0.01 m、寬0.001 m 和高0.04 m.模擬過程采用FVM-DEM 耦合模擬流程,連續相采用有限體積法,離散單元法追蹤顆粒并根據牛頓定律計算求解.采用ICEM 軟件對流體域建模并進行網格劃分,幾何模型導入EDEM 軟件,將修正P-E 磁化模型編寫為動態鏈接庫文件導入,并采用DEM 方法計算離散顆粒運動參數;網格導入FLUENT 軟件,采用FVM 方法計算流體域參數,通過自定義函數實現FLUENTEDEM 雙平臺耦合仿真.將EDEM 時間步長設置為1×10-6s,Fluent 時間步長設置為2×10-5s,兩者的耦合時間步長差為20 倍,節省了計算機資源,提高仿真精度,保持流程時間同步,輸出參數時間同步.

2.2 網格無關性驗證

參考顆?？倓幽苓x取網格尺寸,選取的網格尺寸能夠保持系統內較低的總動能,即可認為網格已經符合要求,其中網格無關性檢驗結果見表1 所示.

表1 網格無關性檢驗表Table 1 Grid independence test table

由表1 可知,采用4R 至3R 所得到的模擬結果較為穩定,對于離散顆粒模擬過程,網格的主要作用是用于顆粒的檢索,過小的網格不利于對顆粒的追蹤,易產生較大誤差.基于求解器的計算穩定性、曳力計算的精度等問題,一般要求流體域網格尺寸大于顆粒至少3 倍.因此綜合多方面的因素考慮選取3R 的網格尺寸,網格數量為248,該網格較為精細且同時滿足雙平臺軟件使用的需求.

2.3 模型驗證

修正P-E 磁化模型的驗證,選取Pinto-Espinoza[18]的實驗數據和Ke 等[20]所使用P-E 磁化模型模擬數據作為對比數據.采用FVM-DEM 耦合模擬方法對雙顆粒運動進行模擬驗證,模擬鐵磁性顆粒在不同條件下所受到的磁化力對運動的影響.模擬參數如表2 所示.

表2 模擬參數Table 2 Analog parameters

選取單個鐵磁性顆粒作為研究對象,模擬當施加垂直方向磁場時,另一個鐵磁性顆粒以固定的顆粒間距圍繞其轉動(顆粒對偏角θ),記錄顆粒所產生的磁化力.選取對應3 種工況,研究鐵磁性顆粒所受磁化力,模擬工況如表3 所示.

表3 模擬工況表Table 3 Simulated working condition table

在不同的工況條件下,分別進行模擬并記錄觀測顆粒所受到的磁化力的作用效果,將模擬結果進行統計并得到磁化力變化圖如圖2 所示.

圖2 磁化力變化圖Fig.2 Magnetization force change diagram

由圖2(a) 和圖2(b) 可知,鐵磁性顆粒在B0=0.02 T 的條件下,不同θ偏角時顆粒所受磁化力.修正P-E 模型的磁化力的計算結果與實驗數據幾乎保持一致,偏差保證在5%以內,且比原模型計算得到的磁化力更貼近實驗數據,因此該模型具有可靠性.

為保證鐵磁性顆粒在磁場環境下的充分運動,通過檢測顆?？倓幽芸梢源_定顆粒是否達到穩定狀態,選取工況C1 顆?？倓幽茏兓瘓D如圖3 所示.

圖3 C1 總動能Fig.3 C1 total kinetic energy

由圖3 可知,總動能初始時刻突然增高,隨流動時間而逐漸降低,該過程中伴隨部分尖峰.由于磁場將顆粒磁化,顆粒受到的磁化力大小與方向,和顆粒對的空間位置有著密切的相關,盡管部分顆粒達到穩定狀態,但細微的變化都會導致顆粒的運動變化,因此當系統內總動能降低到最高動能的5%時,即可認為系統已經達到穩定狀態.

對于研究不同磁場方向及不同配比條件下的鐵磁性顆粒與惰性顆粒在磁場中的運動特性,模擬的各項參數見表4.

表4 模擬參數表Table 4 Analog parameter table

鐵磁性顆粒與惰性顆粒的摻混,采用不同組分配比時,會影響顆粒在磁場中的運動狀態,因此引入鐵磁性顆粒配比變量,公式如下

其中,Nm為鐵磁性顆粒數量,NA為注入顆粒數量.選取3 種不同鐵磁性顆粒占比和磁場方向設計模擬工況,選取磁感應強度B0=0.02 T,其中模擬工況見表5.

表5 模擬工況表Table 5 Simulated working condition table

3 結果與分析

3.1 磁場方向對多組分顆粒運動模擬與分析

基于各工況的模擬結果,選取模擬流動時長為1 s,由于磁場中惰性顆粒的存在,系統無法快速達到穩定狀態,提取0,0.2,0.4,0.6,0.8 和1 s 時刻的顆粒分布圖.為了方便區分鐵磁性顆粒與惰性顆粒,鐵磁性顆粒在分布圖中為藍,惰性顆粒為紅.當鐵磁性顆粒與惰性顆粒占比為1:1(K=50%)時,系統內各時刻顆粒分布圖如圖4 所示.

由圖4(a)和圖4(d)可知,多組分顆粒施加垂直或水平方向的磁場,惰性顆粒未形成聚鏈,鐵磁性顆粒形成沿磁場方向的部分聚鏈,垂直磁場方向保持排斥.但由于惰性顆粒與鐵磁性顆粒的碰撞,使得鐵磁性顆粒的成鏈速度大幅下降,并難以形成較長的聚鏈.

由圖4(b)和圖4(c)可知,多組分顆粒施加β=30°或β=60°的磁場,僅鐵磁性顆粒形成沿磁場方向的聚鏈,垂直磁場方向上相互排斥,其中鐵磁性顆粒夾帶著部分惰性顆粒沿著磁場方向發生偏轉,其偏轉角度與磁場偏角相等,鐵磁性顆粒形成含有一定偏角的聚鏈狀態,需要更長的流動時間,且顆粒成鏈的長度有所下降.

可以發現,惰性顆粒并不會在磁場作用下形成聚鏈,但由于鐵磁性顆粒與惰性顆粒的碰撞,使得部分惰性顆粒夾雜于鐵磁性顆粒聚鏈之間的間隙中.得出結論,多組分顆粒在不同磁場中的運動,由于惰性顆粒的存在,使得鐵磁性顆粒沿著磁場的方向成鏈速度有所下降,形成的聚鏈長度有所下降,部分惰性顆粒夾雜在顆粒聚鏈中.

其中選取各工況1 s 時刻,提取鐵磁性顆粒和惰性顆粒的速度矢量,得到系統內顆粒速度矢量圖如圖5 所示.

圖5 1 s 時刻顆粒速度矢量圖Fig.5 Particle velocity vector diagram at 1 s

由圖5 可知,工況C5 鐵磁性顆粒與部分惰性顆粒的速度矢量指向于分層中心;工況C6 和C7 鐵磁性顆粒與部分惰性顆粒,左側速度矢量垂直向上,右側速度矢量垂直向下,局部顆粒速度矢量形成環形,鐵磁性顆粒與部分惰性顆粒沿著順時針方向逐漸形成沿磁場方向的偏轉;工況C8 局部顆粒速度矢量形成環形.得出結論,多組分顆粒系統下,部分惰性顆粒會被鐵磁性顆粒包裹與夾帶,處于鐵磁性顆粒聚鏈的間隙中,其余惰性顆粒游離在聚鏈外.

進一步提取以上4 種工況不同時刻的總動能,得到顆?？倓幽茏兓€如圖6 所示.

圖6 總動能變化曲線Fig.6 Total kinetic energy change curve

由圖6 可知,當K=50%時,4 種工況所具有的初始能量幾乎相等,工況C5 與C8 采用垂直與水平方向磁場,其總動能曲線下降段的斜率,即顆粒成鏈速度相比含磁場傾角的工況C6 與C7 高,并且曲線下降速度更快,系統總動能更早趨近于0,更為迅速地達到穩定狀態.得出結論,當顆粒配比為1:1 時,施加不同方向的磁場,鐵磁性顆粒初始狀態下能量不變,且當施加水平與豎直方向磁場,多組分顆粒系統達到穩定速度最快.

3.2 鐵磁性顆粒配比對顆粒運動模擬與分析

采用3 種不同鐵磁性顆粒與惰性顆粒的配比,分別研究1:2,1:1 和2:1 的顆粒配比(即K分別為33%,50%和66%)下,顆粒在不同磁場中的運動與分布,選取4 組對照工況,提取其1 s 時刻的顆粒分布圖.組合相同磁場偏角的工況,其1 s 時顆粒分布圖如圖7 所示.

圖7 1 s 時刻顆粒分布圖Fig.7 Particle distribution map at 1 s

由圖7(a)可知,工況C1,C5 和C9 采用垂直磁場,鐵磁性顆粒占比逐漸升高.各工況均出現了顆粒聚鏈現象,C1 形成較短聚鏈,C5 聚鏈數量增多且開始出現分層現象,C9 形成了較長聚鏈且出現明顯分層現象.

由圖7(b)和圖7(c)可知,工況C2,C6 和C10 采用β=30°傾角磁場,工況C3,C7 和C11 采用β=60°傾角磁場,隨著鐵磁性顆粒占比逐漸升高,鐵磁性顆粒沿磁場方向的聚鏈增多,并形成長鏈,而且鐵磁性顆粒對惰性顆粒的束縛能力逐漸增強,更為緊密地將惰性顆粒包裹在鐵磁性顆粒的聚鏈中.

由圖7(d)可知,工況C4,C8 和C12 采用水平磁場.分析不同配比下鐵磁性顆粒的聚鏈變化,得出同樣結論,隨著鐵磁性顆粒占比增加,聚鏈數量與成鏈長度同樣有所增加,鐵磁性顆粒對惰性顆粒的束縛能力逐漸增強,將惰性顆粒更緊密地包裹在顆粒聚鏈中.

提取以上12 種工況不同時刻的總動能,得到顆?？倓幽茏兓€如圖8 所示.

圖8 總動能變化曲線Fig.8 Total kinetic energy change curve

由圖8 對各組內部數據進行對比,能夠直觀得到結論,當鐵磁性顆粒占比越高,顆粒初始狀態所具有的能量均有所增加.觀察曲線的變化,通過總動能下降過程中的斜率,進而判斷成鏈速度達到穩定快慢,得到結論,鐵磁性顆粒占比越高,多組分顆粒系統越容易達到穩定,顆粒的成鏈速度也就越快.且由圖8(b)和圖8(c)可知,含有磁場偏角的多組分顆粒,總動能變化劇烈,運動中較難達到穩定狀態.

認為總動能降低為最高動能5% 則達到穩定,經過整理得到穩定時刻直方圖如圖9 所示.

圖9 各工況穩定時刻直方圖Fig.9 Histogram of stable time of each working condition

圖9 中工況C1,C2,C3,C4,C6,C7 和C9 未在1 s 內達到穩定狀態,系統仍具有波動.對比圖中工況C4,C8 和C12 數據可知,當采用水平方向磁場時,鐵磁性顆粒占比越大,穩定所需時間從1 s 以上縮短至0.2 s,顆粒達到穩定狀態的速度顯著增快.

4 結論

本文采用FVM-DEM 耦合模擬方法,結合修正P-E 模型,模擬研究鐵磁性顆粒在惰性顆粒摻混的條件下,不同磁場方向與不同組分配比對鐵磁性顆粒運動的影響,通過模擬得到各個工況的數據參數,整理得到顆粒分布圖、顆粒速度矢量圖和顆?？偰芰孔兓瘓D,對模擬結果進行分析,總結了多組分顆粒在不同工況下的運動規律.為鐵磁性顆粒在不同磁場中的運動,提供了理論模型與模擬參考,得出主要結論如下.

(1) 在多組分顆粒系統中,施加不同傾角的磁場時,鐵磁性顆粒仍會沿磁場方向形成聚鏈,但其顆粒聚鏈長度與成鏈速度相比垂直磁場有所下降,且部分惰性顆粒會被鐵磁性顆粒包裹與夾帶,處于鐵磁性顆粒聚鏈的間隙中;鐵磁性顆粒占比為1:1 時,鐵磁性顆粒具有的初始能量不受磁場方向的影響,且當施加水平與豎直方向磁場時,多組分顆粒系統達到穩定速度最快.

(2) 在多組分顆粒系統中,隨著鐵磁性顆粒占比升高,鐵磁性顆粒具有的初始能量將有所升高、聚鏈數量與成鏈長度將有所增加、鐵磁性顆粒對惰性顆粒的束縛能力逐漸增強;當施加含有傾角的磁場時,隨著鐵磁性顆粒占比升高,鐵磁性顆粒達到穩定狀態需要的時間逐漸縮短.

(3) 在多組分顆粒系統中,施加水平磁場時,可以通過增大鐵磁性顆粒占比有效提升穩定速度,使系統更快趨于穩定;但對于施加磁場偏角的工況,較難通過改變鐵磁性顆粒占比縮短穩定所需時間.