磁性液體兩相界面演變特性的數(shù)值研究

施東曉,畢勤成,周榮啟

(西安交通大學動力工程多相流國家重點實驗室, 710049, 西安)

磁性液體兩相界面演變特性的數(shù)值研究

施東曉,畢勤成,周榮啟

(西安交通大學動力工程多相流國家重點實驗室, 710049, 西安)

為了克服VOF,Level Set等方法在磁性液體兩相流動研究中存在的不足,采用界面追蹤法(VOSET, coupled Volume-of-Fluid and Level Set)捕捉相界面,將磁場力作為源項添加到流體運動的動量方程中,建立了磁性液體流場與磁場耦合的數(shù)值方法。推導了麥克斯韋方程的二維理論解并用于驗證數(shù)值解,兩者吻合良好,驗證了磁場求解的正確性;同時通過模擬磁性液體液滴的平衡形狀驗證了本數(shù)值方法的準確性;進而分別研究了均勻磁場作用下磁性液體液滴振蕩及非磁性氣泡在磁性液體中上升的問題。研究結果表明,液滴或氣泡均會沿磁場方向拉長,增大磁場強度或磁化率將導致更顯著的變形;磁感線在相界面出現(xiàn)彎曲,彎曲方向指向磁導率較大的介質。

磁性液體;兩相流;數(shù)值模擬;界面追蹤法

磁性液體(ferrofluids或magnetic fluids)是含有顆粒直徑10nm左右的單磁疇磁性顆粒穩(wěn)定膠體懸浮液[1]。與磁流變流體在外加強磁場中會發(fā)生固化不同,包覆在顆粒表面的活性劑分子使得磁性液體仍可以保持良好的流動性。憑借獨特的磁性與流動性,磁性液體在真空密封、揚聲器冷卻、生物醫(yī)學工程等領域已有了一定的實際應用[2]。在這些應用中,往往涉及到了帶有相界面的磁性液體兩相流動。由于外部磁場施加在相界面上麥克斯韋應力的作用,磁性液體的自由表面流動與普通流體存在很大的差異。研究表明,磁性液體液滴或磁性液體中的非磁性氣泡均會沿磁場方向拉長[3-4],施加一個非均勻磁場可以控制磁性液體液滴的運動軌跡[5-6],這使得遠程非接觸操控液滴成為可能。建立準確的數(shù)理模型有助于進一步揭示磁場和流場的耦合機制,拓展磁性液體的應用。

氣泡或液滴的運動是兩相流中的基本問題,目前已有學者針對磁性液體氣泡/液滴的行為開展了相關的研究。Korlie等采用VOF方法模擬了氣泡在磁性液體中的上升過程[7];Zhu等運用Level Set方法研究了磁性液體液滴在超疏水表面上的非線性變形[8];Liu等采用Particle Level Set方法模擬了微流體流動聚焦結構中磁性液體液滴的生成過程[9]。VOF和Level Set方法因為可以處理帶有復雜界面拓撲變化的兩相流動問題而被廣泛采用,但這兩種方法均存在各自的缺點。VOF計算界面曲率及法向不夠準確,且求得的界面附近的物性分布是不連續(xù)的,將影響表面張力和磁力的計算精度;Level Set則存在質量損失問題,難以模擬尖銳的界面。Particle Level Set雖然能夠在一定程度上解決質量損失問題,但該方法執(zhí)行較為復雜,通用性不強。界面追蹤法(VOSET)[10]是近幾年剛提出的一種耦合VOF和Level Set方法的新型界面捕捉技術,不僅繼承了VOF方法的高質量守恒性,還結合了Level Set函數(shù)光滑性好的優(yōu)點,可以提高界面曲率的計算精度并光順界面附近突變的物理量。VOSET只需求解流體體積函數(shù)的輸運方程,而Level Set函數(shù)則是通過幾何方法確定,極大地簡化了計算過程。

本文首次將VOSET方法與磁場力模型相結合,建立了模擬均勻磁場中磁性液體不可壓縮兩相流動的二維數(shù)值模型?？刂品匠贪o磁場非導電流體的麥克斯韋方程,以及耦合表面張力和磁力的動量方程。為了驗證麥克斯韋方程求解的準確性,首先推導了適用于二維坐標系下的磁勢拉普拉斯方程的理論解,并將數(shù)值解與其進行對比。同時通過模擬磁性液體液滴的平衡形狀驗證了本數(shù)值方法的正確性。隨后根據(jù)本文建立的數(shù)值模型分別研究了均勻磁場作用下磁性液體液滴振蕩及氣泡在磁性液體中上升的問題。

1 數(shù)學模型與數(shù)值方法

1.1 控制方程

假設納米磁性顆粒在基液中為均勻的單分散狀態(tài),且相互之間不發(fā)生作用,在對磁性液體宏觀行為的研究中,可以將磁性液體當作性質單一的均勻流體進行處理[7,9]。磁場對磁性液體界面變形的影響來源于作用在相界面上的麥克斯韋應力,該應力可以簡化為界面上的法向應力。考慮不可壓縮兩相流體的流動,表面張力和磁力存在于相界面,連續(xù)方程及動量方程可描述為

(1)

(2)

式中:u、ρ、η和p分別為流體的速度、密度、黏度和壓力;Fσ為表面張力。方程(2)右邊的最后一項為磁力項,采用與表面張力相同的處理方式,將磁力以源項的形式添加到動量方程中。麥克斯韋應力張量τm有如下形式

(3)

式中:μ為磁導率;H為外加磁場;I為單位張量。

1.2 界面追蹤法

VOSET采用VOF方法捕捉相界面,保持了兩相間的質量守恒;構造界面附近的Level Set符號距離函數(shù)φ用來計算界面曲率、法向和光滑不連續(xù)的物性,克服了VOF方法計算曲率不精確和光滑物性差的缺點。

在VOF方法中,采用流體體積函數(shù)C表示計算單元中流體占據(jù)單元空間的體積分數(shù),定義C=0表示網(wǎng)格內全為連續(xù)相,C=1表示網(wǎng)格內全為離散相,0

(4)

方程(4)的求解和相界面的重構采用了Youngs的PLIC算法[11]。

根據(jù)求得的流體體積函數(shù),通過迭代幾何計算方法求出界面兩側各3倍網(wǎng)格寬度區(qū)域內的Level Set函數(shù),步驟如下:①在整個計算區(qū)域內為Level Set函數(shù)賦初值;②標識相界面附近的計算單元;③計算標識區(qū)域內單元(i,j)到界面的最短距離d,根據(jù)下式確定符號距離函數(shù)的值

(5)

④用求得的Level Set函數(shù)φ重新求界面法向n=φ/|φ|,根據(jù)新的法向再次構造界面,重復步驟③、④,詳細的計算過程參考文獻[10]。

采用Level Set函數(shù)來光滑密度、黏度和磁導率

ρ(φ)=ρd+(ρc-ρd)Hε(φ)

(6)

η(φ)=ηd+(ηc-ηd)Hε(φ)

(7)

μ(φ)=μd+(μc-μd)Hε(φ)

(8)

式中:下標d和c分別表示離散相和連續(xù)相;Hε(φ)為光滑Heaviside函數(shù)[12],其表達式為

(9)

本文取系數(shù)ε=1.5Δ,其中Δ為網(wǎng)格寬度。

表面張力由連續(xù)表面力(CSF)模型給出,寫成Level Set函數(shù)形式[12]有

Fσ=-σκ(φ)Hε(φ)

(10)

式中:σ為表面張力系數(shù);κ為界面曲率,由下式得出

κ(φ)=·

(11)

1.3 磁場力模型

一般認為磁性液體為非導電流體,且流體中不存在位移電流,描述磁場的麥克斯韋(Maxwell)方程組則可簡化為

(12)

施加外磁場時,磁性液體內部的磁性顆粒磁矩沿外磁場方向排列,使得磁性液體對外顯磁性,定義單位體積內顆粒磁矩的矢量和為磁化強度M。磁感應強度B、磁場強度H以及磁化強度M之間滿足本構關系

B=μ0(H+M)

(13)

(14)

引入標量磁勢ψ,定義H=-ψ,將方程(12)的求解轉化為求解關于磁勢的拉普拉斯方程·(μψ)=0

(15)

在外加磁場作用下,流體中產生的磁界面力可以由麥克斯韋應力張量給出Fm=·τm,采用與表面張力相同的處理方式,將表面力轉化為體積力,那么單位體積磁力的表達式為[1]

Fm=-μ0H,TdH+μ0MH

(16)

磁化強度M(μ,H,T)為外磁場和溫度的函數(shù),將式(16)展開得

Fm=-μ0MdH+μ0υH,TdH+

μ0MH

(17)

式(17)右邊第一項表示流體-磁壓力項,即

(18)

(19)

此外,式(17)右邊第二項為磁致伸縮壓力項,即

(20)

將M=(μ-μ0)H/μ0及比容υ=1/ρ代入上式,得ps=0。

基于以上分析,磁力公式(16)最終簡化為

(21)

1.4 數(shù)值方法

在二維均勻直角網(wǎng)格上采用有限容積法離散控制方程,磁勢方程(15)的離散采用二階精度的中心差分法,方程(1)和(2)的求解在交錯網(wǎng)格上執(zhí)行,壓力與速度的耦合采用SIMPLER算法。一個時間步長內的計算過程如下:①由t時刻的流體體積函數(shù)Cn根據(jù)VOSET方法求出Level Set函數(shù)φn;②求解方程(15)獲得磁場分布Hn;③ 由求得的φn和Hn計算表面張力和磁力,代入方程(2),求出t+1時刻的速度場un+1;④根據(jù)新的速度場un+1,采用Youngs的PLIC界面重構技術構造下一時刻的界面Cn+1。

2 結果與分析

2.1Maxwell方程求解與驗證

導磁圓球在均勻磁場中的磁勢方程存在理論解[1],可用于驗證Maxwell方程的數(shù)值解。球坐標系下磁勢理論解表征的是三維磁場分布,針對二維情況,圓球模型不再適用,應采用無限長圓柱體模型,如圖1所示,計算區(qū)域為LX×LY,中心圓柱體半徑為R,內部磁導率為μ1,外部磁導率為μ2,施加豎直向上的均勻磁場H0。

圖1 計算模型示意圖

(22)

式中:r為空間中任意一點到圓心的距離。方程(22)的解為

(23)

根據(jù)磁勢定義可求得磁場分布

(24)

Maxwell方程的邊界條件[1]為

n·(B1-B2)=0

n×(H1-H2)=0

(25)

式中:n為界面法向;B1、B2、H1和H2分別表示界面兩側的磁感應強度和磁場強度矢量。由式(25)可知B場在界面法向連續(xù),即B1n=B2n;H場在界面切向連續(xù),即H1t=H2t。當r→∞時,H=H0j(j=ersinθ+eθcosθ),因此E=-H0。根據(jù)界面邊界條件,求得系數(shù)A和D

(26)

采用有限容積法離散磁勢方程,對于二維直角坐標,磁勢方程(15)可寫成

(27)

選取計算參數(shù):LX=0.02m,LY=0.02m,R=0.002m,μ2=μ0,μ1=2μ0,H0=1 kA·m-1,網(wǎng)格200×200。圖2對比了y軸上的磁場值,理論解求得圓內磁場值為0.667 kA·m-1,數(shù)值解為0.662kA·m-1,兩者基本吻合,表明Maxwell方程求解是正確的。

圖2 沿y軸(x=0.01 m)的磁場強度

圖3給出了圓柱周圍磁場及磁感線的分布,圓柱內部磁場分布均勻,磁感線與外磁場同向;圓柱內外磁導率的突變導致磁感線在界面附近發(fā)生扭曲,而在遠離圓柱的區(qū)域又恢復為均勻。

圖3 磁場及磁感線分布

2.2 磁性液體液滴的平衡形狀

為驗證本文建立的模型及方法,首先計算了初始靜止的磁性液體液滴在均勻磁場中的平衡形狀問題,并將計算結果與Flament等人的實驗結果[13]進行對比。半徑R=0.001 m的液滴位于8R×8R的區(qū)域中心,磁化率為2.2,表面張力為3.07 mN·m-1,不計重力,壁面采用無滑移速度邊界條件,計算網(wǎng)格取128×128。圖4給出了不同大小磁場作用下液滴的平衡形狀及高寬比與實驗結果的對比。定義液滴在平行和垂直于磁場方向上的直徑的比值E來表征液滴的變形程度[14]。圖4中工況1～4所對應的磁場大小分別為1.2、2.4、2.9及3.7 kA·m-1?？梢钥闯?計算結果與實驗結果吻合較好,證明本文的數(shù)值算法可有效地計算磁性液體的兩相流動問題。

圖4 計算結果與實驗結果[13]的對比

2.3 磁性液體液滴振蕩

下面模擬磁性液體液滴在均勻磁場中的振蕩過程。計算區(qū)域為0.08 m×0.08 m,半徑為0.01 m的圓形液滴位于區(qū)域中心,液滴及周圍空氣的物性參數(shù)為ρd=103kg·m-3,ρc=1 kg·m-3,ηd=10-3Pa·s,ηc=10-5Pa·s,μd=10μ0,μc=μ0,σ=0.1 N·m-1。外加豎直向上的均勻磁場H0=2kA·m-1,不考慮重力,四周設為無滑移固壁,計算網(wǎng)格取112×112。

施加磁場后,麥克斯韋應力在液滴沿磁場方向的兩極取最大[7],為保持法向應力平衡,液滴沿磁場方向拉伸成扁長形狀,而表面張力阻礙了液滴的進一步拉伸,因此在磁力和表面張力的共同作用下液滴沿豎直方向振蕩變形,如圖5a～5f所示。圖5g～5i為對應于圖5a～5c時刻的磁感線分布,可以看出磁感線隨液滴形狀的變化而改變,且磁感線朝著液滴方向彎曲,這是由于液滴磁導率大于空氣磁導率。

(a)0s (b)0.07 s (c)0.16 s

(d)0.24s (e)0.30s (f)0.37 s

(g)0s (h)0.07 s (i)0.16 s

2.4 磁性液體中氣泡上升

考慮8R×16R的計算區(qū)域,半徑R=0.004m的氣泡初始位于(4R,4R),均勻磁場豎直向上,網(wǎng)格取112×224。采用Eo數(shù)和Mo數(shù)描述氣泡的運動特性

(28)

(a)無磁場 (b)豎直磁場(Bom=4)

圖7為0.16 s時不同磁場中的氣泡形狀及磁感線分布,Bom=2時磁場較弱,氣泡仍為扁平的橢圓形;隨著磁場強度的增大,氣泡逐漸沿磁場方向變形,且磁場越強變形越明顯。磁感線在氣泡附近發(fā)生扭曲,并背離氣泡方向彎曲,與磁性液體液滴的情況相反。

(a)Bom=2 (b)Bom=4

(c)Bom=6 (d)Bom=8

另外,由磁力公式(21)可知,磁力的大小與磁化率成正比,因此氣泡的運動行為還受磁化率的影響。圖8展示了相同磁場作用下不同磁化率時的氣泡形狀(Bom=4,t=0.16 s),可以看出氣泡的變形程度隨著磁化率的增加而增大。

(a=1 (b=2

(c=3 (d=4

圖9給出了氣泡的高寬比、最終速度與Bom及磁化率的關系。由于增大磁場強度(Bom)或磁化率均會增大磁力,從而促進了氣泡沿磁場方向上的變形。氣泡的高寬比越大,受到的阻力越小,上升得越快,因此氣泡的形狀影響著氣泡的上升速度。

(a)氣泡的高寬比及速度與Bom的關系

(b)氣泡的高寬比及速度的關系

3 結 論

本文對磁場力模型進行了詳細推導,采用與表面張力相同的處理方法,給出了磁力的最終表達式,并作為源項添加到動量方程中,結合VOSET方法,建立了模擬磁性液體不可壓縮兩相流動的二維數(shù)值模型。推導了Maxwell方程的二維理論解,并分別計算了均勻磁場中磁性液體液滴振蕩與氣泡在磁性液體中上升的問題,得出以下結論。

(1)Maxwell方程的數(shù)值解與推導的理論解吻合較好,表明本數(shù)值方法求解Maxwell方程的正確性。

(2)麥克斯韋應力在磁性液體液滴沿磁場方向的兩極取最大,為保持法向應力平衡,液滴被拉長為扁長的形狀,而表面張力阻礙了液滴的進一步拉伸,液滴最終在磁力和表面張力的作用下振蕩變形。

(3)磁性液體中的非磁性氣泡在磁場作用下克服表面張力沿磁場方向拉伸變形,增大磁場強度或磁導率均會增大氣泡的高寬比,減小了氣泡在運動方向上的阻力,提高了上升速度。

(4)磁感線隨著液滴或氣泡形狀的變化而變化,并在界面附近(磁導率突變的位置)發(fā)生扭曲,彎曲方向指向磁導率較大的介質。

(5)本文建立的磁性液體流場與磁場耦合的數(shù)值方法,可以為帶有相界面的磁性液體流動研究提供理論基礎。

[1] ROSENSWEIG R E. Ferrohydrodynamics [M]. Cambridge, UK: Cambridge University Press, 1985.

[2] RINALDI C, CHAVES A, ELBORAI S, et al. Magnetic fluid rheology and flows [J]. Current Opinion in Colloid & Interface Science, 2005, 10(3/4): 141-157.

[3] ARKHIPENKO V I, BARKOV Y D, BASHTOVOI V G. Shape of a drop of magnetized fluid in a homogeneous magnetic field [J]. Magnitnaia Gidrodinamika, 1979, 14: 131-134.

[4] ISHIMOTO J, OKUBO M, KAMIYAMA S, et al. Bubble behavior in magnetic fluid under a nonuniform magnetic field [J]. JSME International Journal: Series B Fluids and Thermal Engineering, 1995, 38(3):

382-387.

[5] NGUYEN N T, NG K M, HUANG X. Manipulation of ferrofluid droplets using planar coils [J]. Applied Physics Letters, 2006, 89(5): 052509.

[6] PROBST R, LIN J, KOMAEE A, et al. Planar steering of a single ferrofluid drop by optimal minimum power dynamic feedback control of four electromagnets at a distance [J]. Journal of Magnetism and Magnetic Materials, 2010, 323(7): 885-896.

[7] KORLIE M S, MUKHERJEE A, NITA B G, et al. Modeling bubbles and droplets in magnetic fluids [J]. Journal of Physics: Condensed Matter, 2008, 20: 204143.

[8] ZHU G-P, NGUYEN N T, RAMANUJAN R V, et al. Nonlinear deformation of a ferrofluid droplet in a uniform magnetic field [J]. Langmuir, 2011, 27(24): 14834-14841.

[9] LIU J, YAP Y F, NGUYEN N T. Numerical study of the formation process of ferrofluid droplets [J]. Physics of Fluids, 2011, 23: 072008.

[10]SUN D L, TAO W Q. A coupled volume-of-fluid and level set (VOSET) method for computing incompressible two-phase flows [J]. International Journal of Heat and Mass Transfer, 2010, 53(4): 645-655.

[11]YOUNGS D L. Time-dependent multi-material flow with large fluid distortion [M]∥MORTON K W, BIANES M J. Numerical Methods for Fluid Dynamics. New York, USA: Academic Press, 1982: 273-285.

[12]SUSSMAN M, FATEMI E, SMEREKA P, et al. An improved level set method for incompressible two-phase flows [J]. Computers & Fluids, 1998, 27(5/6): 663-680.

[13]FLAMENT C, LACIS S, BACRI J C, et al. Measurements of ferrofluid surface tension in confined geometry [J]. Physical Review: E, 1996, 53: 4801.

[14]劉艷艷, 李彥鵬, 朱婷婷. 表面活性劑對中尺度氣泡形狀及速度的調控研究 [J]. 西安交通大學學報, 2011, 45(10): 93-97.

LIU Yanyan, LI Yanpeng, ZHU Tingting. Study on modulating shape and velocity of meso-scale bubble using surfactants [J]. Journal of Xi’an Jiaotong University, 2011, 45(10): 93-97.

(編輯 荊樹蓉)

NumericalInvestigationontheEvolutionCharacteristicsofTwo-PhaseInterfaceinFerrofluids

SHI Dongxiao,BI Qincheng,ZHOU Rongqi

(State Key Laboratory of Multiphase Flow in Power Engineering, Xi’an Jiaotong University, Xi’an 710049, China)

To overcome the disadvantages of VOF and Level Set methods in the research of ferrofluid two-phase flows, a numerical method coupling the flow field with magnetic field of ferrofluids was developed, where the evolution of interface was captured by the VOSET (coupled Volume-of-Fluid and Level Set) method, and the magnetic force was added into the momentum equation as a source term. To verify the accuracy of solving magnetic fields, a 2D analytical solution of Maxwell equations was derived and used to validate the numerical solution. The numerical results were in good agreement with the theoretical values. Then, the equilibrium shape of a ferrofluid droplet was simulated which verified the validity of our numerical method. Based on the proposed method, the droplet oscillation and rising bubble problems in the presence of magnetic fields were numerically investigated. The results showed that both the droplet and bubble were stretched in the direction of magnetic field, and the increase in either magnetic field strength or magnetic susceptibility would lead to a more significant deformation. Note that the magnetic field lines were distorted in the vicinity of interface with the bending direction pointing to the medium with a higher susceptibility.

ferrofluids; two-phase flow; numerical simulation; interface tracking method

2014-02-04。

施東曉(1987—),男,博士生;畢勤成(通信作者),男,教授,博士生導師。

國家自然科學基金資助項目(11072189);國家“863計劃”資助項目(2008AA05Z417)。

時間:2014-06-13

10.7652/xjtuxb201409021

O359.1

:A

:0253-987X(2014)09-0123-07

網(wǎng)絡出版地址:http:∥www.cnki.net/kcms/detail/61.1069.T.20140613.1457.005.html