基于WENO格式有限體積法的鐵磁流體兩相流相場方法

中圖分類號： 0359÷.1 文獻標志碼：A 文章編號：1000-582X（2025）04-067-17

A finite volume-based phase field method for two-phase ferrofluid flows

ZHANG Shaosong'， ZHANG Liangqi'， CHENG Liming'， WANG Xiaoshuang2， XIAO Yao'， ZENG Zhong' （1.College of Aerospace Engineering，Chongqing University， Chongqing 40o044， P. R. China； 2.ShanghaiChongqing Institute of Artificial Intelligence， Chongqing 401332， P. R. China）

Abstract：This paper presents a finite volume method based on the weighted esentially non-oscillatory（WENO） scheme to develop a phase field method for simulating two-phase ferrofluid flows.The incompressible Navier

Stokes quations is used to describe fluide flow，the Cahn-Hillard equation is adopted to capture interfacial motion of two-phase flow， and the Maxwellequation is used to describe external magnetic field distribution.At the same time， adding Kelvin force and surface tension to the fluid flow control equation to achieve the description of interface dynamic behavior by magnetic field.To address the chalenges posed bythe fourth-order nonlinear diffusion terms，the Cahn-Hilliard equation is decomposed into two Helmholtz equations.The fifthorder WENO scheme is employed to handle the convection term， enhancing computational accuracy and mitigating numerical oscillations. Validation through Zalesak's disk problem shows that the proposed method achieves higher phase interface capture accuracy compared to existing references，while maintaining performance comparable to high-precision phase field methods.The method isapplied to investigate droplet shear deformation， revaling its capability to capture more satellte droplets. Moreover，research on the shear deformation of ferromagnetic fluid droplets under the influence of magnetic fields and with lower capillary numbers indicatesthat the magnetic interfacial force favors droplet deformation when the external magnetic field direction aligns with the hydrodynamic deformation.Furthermore， increasing the magnetic field intensity leads to droplet splitting. Conversely， when the magnetic field is nearly perpendicular to the deformation direction， a low-intensity field alters the deformation trajectory，while a high intensity magnetic field enforces deformation along the magnetic field direction.

Keywords： phase-field method； finite volume method； ferrofluid droplets； magnetic field control

鐵磁流體憑借其優越的變形能力和磁場可控性而受到廣泛應用。最具代表性的應用是鐵磁流體液滴可用作反應平臺或載體，在芯片實驗室中實現高效可控的樣品輸運。此外，鐵磁流體液滴在外加磁場的精確控制下能分裂成多個子液滴，再按需重新融合，被視為能實現復雜功能的微型磁性軟機器人，在醫療保健、生物工程和微操作等領域展現了巨大應用潛力，吸引了大量研究者的興趣[。

為了深入研究鐵磁流體的界面動力學行為，已經發展出多種多相鐵磁流體流動數值方法。Gao等采用流體體積（VOF）方法模擬了微重力真空條件下金屬液滴在磁場區域內的運動情況。Korlie等提出了一種VOF方法來模擬鐵磁流體中氣泡和液滴的運動。CunhaLucas等利用LevelSet方法對在較高雷諾數（ R e= 1）下，磁場誘導的剪切流中鐵磁液滴破碎的現象進行數值模擬研究。研究結果表明，當磁場作用方向與流動方向平行時，會延遲破裂過程，減小衛星液滴的尺寸；相反，如果在垂直于流動方向上施加磁場，則可以通過調整磁場強度來誘導、延遲甚至阻止液滴的破裂。近期有許多學者開始使用相場法來描述多相鐵磁流體流動行為。與前述的2種界面描述方法相比，相場法避免了額外的界面幾何重構和重新初始化過程，而且表面張力計算簡單，并且能夠確保質量和能量的守恒。Nochetto等[開發出含鐵磁流體兩相流的相場模型，能捕獲鐵磁流體的基本流動現象，如Rosensweig不穩定性和鐵磁流體刺猬。Hu等將相場方法應用于鐵磁流體流動，采用守恒的Allen-Cahn方程捕捉界面。他們開發了一種多松弛格子玻爾茲曼方法（MRT-LBM）來求解控制方程。Khan等用新近發展的簡化多相格子玻爾茲曼方法（SMLBM）捕捉鐵磁流體的流體動力學行為，通過數值模擬研究，探討了在均勻磁場作用下鐵磁流體中鐵磁液滴變形和磁流體氣泡融合的動力學過程。

鐵磁流體液滴在微流控技術和生物醫藥領域的前沿應用需要精確控制其界面動力學行為，而對鐵磁液滴運動和變形行為進行數值模擬研究能為理解控制機理和設計控制方式提供必要的理論支持。筆者擬選擇鐵磁流體相場模型為理論模型，應用基于五階WENO格式的有限體積方法統一離散求解控制方程組，以準確描述鐵磁流體液滴的界面動力學行為。Tang等指出，基于Cahn-Hilliard方程相場模型的數值求解存在很大難度，如處理高階非線性擴散項、保持長時間數值計算的精度，以及能量穩定性等，都是需要克服的挑戰。Yue等明確指出相變量的數值擾動會導致液滴體積的非物理收縮，并破壞質量守恒。這里，借鑒了Dong等的工作思路，將四階Cahn-Hilliard方程拆分為2個Helmholtz形式方程，這樣可以克服由四階非線性項帶來的困難。為了提升魯棒性，采用了五階WENO格式來統一處理控制方程的對流項，以提高數值計算的精度，同時避免數值振蕩。此外，參考Huang等的研究，引入了一致、守恒的質量通量來描述相變量擴散對質量和動量輸運的影響，保證相界面處動量輸運的守恒性。采用增量壓力修正方案，實現不可壓縮流的流動速度和壓力的解耦，從而在每一個時間步，只需要求解一系列解耦的速度和壓力離散方程。本文提出的方法采用了統一的控制體構造方式，所有變量都存儲在控制體中心，也就是采用同位網格。采用了Rhie-Chow插值修正面心速度，以抑制壓力振蕩。在采用Zalesak's圓盤問題驗證所提出方法在捕捉界面精度上的優越性后，將其應用于研究在磁場作用下鐵磁流體液滴剪切變形的動力學。

1理論模型與數值方法

1.1 兩相鐵磁流體相場模型

研究對象是包含鐵磁流體的兩相流，其中一相為鐵磁流體，另一相為非磁性流體。兩相鐵磁流體流動的相場模型歸納如下。

1.1.1 相場方程

1.1.3 磁場方程

由于關注重點在鐵磁流體的界面動力學，故忽略了鐵磁流體的內旋效應，而采用非導電鐵磁流體的Maxwell方程\"來描述磁場。

式中： B 和 H 分別為磁感應強度和磁場強度。 B 可以表示為

1.2控制方程的離散

采用同位網格進行離散，即采用統一的控制體構造方式，而且所有變量都存儲在控制中心，對流項采用五階WENO顯式處理。

1.1.1 時間離散

1.1.2 空間算子離散

基于中心差分格式，對控制方程的梯度、散度、拉普拉斯算子進行離散，梯度算子離散為

1.1.3 Cahn-Hilliard方程的離散

將相場方程分解為2個二階Helmholtz方程離散求解[15]

1.1.4Navier-Stokes方程的離散

采用增量壓力修正投影法對Navier-Stokes方程進行解耦求解，將動量方程（5）拆分為方程（31）和方程（32）。

1）預測步。將方程（31）的部分黏性項移到方程左側進行隱式處理，對流項采用五階WENO顯式處理，如下所示。

2）壓力修正步。在方程（32）兩端同時取散度，結合連續性方程（6）推導出增量壓力泊松方程

3）速度插值。采用同位網格進行離散，即采用統一的控制體構造方式，而且所有變量都存儲在控制中心，離散方程（34）時引人Rhie-Chow插值，以避免壓力震蕩。將式（33）在如圖1所示的控制體積 P 中進行離散，得到的各項系數為

2數值驗證

2.1界面捕捉精度的數值驗證

利用經典的Zalesak's圓盤旋轉案例來驗證本文方法在界面捕捉精度方面的效果。理論上，經過一個或者多個周期的旋轉后，圓盤將回到初始位置，因此，通過比較經過周期旋轉后的界面形態與初始界面形態即可反映數值計算界面捕捉精度。算例設置如圖2所示。一個帶凹槽的圓盤放在大小為 的周期區域的中心，取 ，圓盤的半徑和凹槽的寬度分別為80和16，擴散界面區域的厚度設置為 ε=2.0 。圓盤在如下速度場的驅動下做周期性旋轉。

圖4展示了 P e=800 時，不同網格分辨率 N×N 的誤差變化曲線，橫坐標為網格分辨率 N（N=50，100，200 256），縱坐標為誤差 。結果表明，誤差隨網格分辨率呈線性收斂，收斂階數約為2.5。

2.2液滴在剪切流動下的變形

單個液滴在剪切流動中的變形是檢驗兩相流數值方法界面捕捉精度的另一個經典算例，該算例設置如圖5所示。一個半徑為 R 的液滴最初放置在 8R×8R 的計算域中心，網格分辨率為 200×200 。計算域的上下固壁分別以大小為 U 的速度向相反方向移動，左右邊界為周期性邊界條件。為了量化液滴變形，定義液滴的變形參數 D=（L-B）/（L+B），L 和 B 分別為變形液滴在平衡狀態下的長軸和短軸的長度。

圖6展示了不同毛細數下（ "C a=0.1，0.25，0.4，0.5） 平衡態液滴界面形態，可以看出隨著 C a 數的不斷增大，黏性力作用增強，液滴變形程度越大。圖7給出了當前數值計算得到的變形參數隨 C a 的變化，當前方法所得數值結果與Cox變形理論和參考文獻[30]非常吻合。

為驗證當前數值方法捕捉液滴破碎行為的能力，探討了在更高雷諾數 （R e） 條件下的液滴剪切變形。取R e=1 ， C a=0.5 λ=1.2 ， ；半徑為 R 的液滴放置在 16R×8R 的計算域的中心，網格分辨率為200×100 0

圖8展示了不同時刻液滴界面形狀的演化過程，觀察到液滴破裂前的界面變形結果與文獻[9]和文獻[30]的報道是一致的。然而，在液滴破碎以后，本文的相場方法預測出現了4個衛星液滴，而相同條件下，文獻[9]報道只出現了1個衛星液滴，文獻[30]的結果出現了2個衛星液滴。 Hu 等[和Li等[32指出，改善守恒性、減少數值耗散能夠顯著提高相場方法捕捉微小界面變形細節的能力，他們改良的相場方法分別在液滴剪切問題和液滴撞擊液膜問題中捕捉到了更多的衛星液滴。此外，相場方法界面捕捉精度與界面厚度和空間網格分辨率相關。這里首先對當前的數值計算進行 C n 數和網格分辨率的收斂性驗證。如圖9所示，較大 C n 數對應過寬的界面厚度，使微小界面結構變化被抹平，所以在 C n=0.01 時，沒有觀測到衛星液滴；降低 C n 數，數值計算捕捉細小界面變形的能力顯著改善，當 C n=0.008 時出現4個衛星液滴，再進一步降低界面厚度，不再出現更多衛星液滴，所以當前結果是關于界面厚度參數收斂的。此外，如圖10所示，固定 C n=0.0075 ，網格分辨率增大，液滴破碎后仍觀察到4個衛星液滴，所以網格分辨率為 200×100 時數值結果已經實現收斂。

2.3靜態磁場作用下液滴的變形

為驗證當前方法描述磁場作用的準確性，模擬了不同磁場強度條件下的鐵磁流體液滴變形。在這種情況下，將半徑為1的鐵磁流體液滴放置在 16×16 的矩形區域的中心。鐵磁流體和非磁性環境流體的材料性質與Flament實驗研究類似。當施加垂直向上方向的均勻磁場時，開爾文力會在2種流體界面處由于磁場作用產生，并打破表面張力與壓力梯度的平衡，鐵磁流體液滴發生沿磁場方向的變形。選用的參數為黏度比 λ= 1，磁導率比 ，網格分辨率為 200×200 0

圖11（a）描述了在不同磁場強度 H（H=1.2，2.4，2.9，3.7，5.5kA/m） 下，鐵磁流體液滴的平衡態界面形狀。當施加豎直方向的均勻磁場時，液滴會在磁界面力的作用下沿磁場方向伸長，而表面張力的作用則努力使液滴恢復到圓形。當液滴開始變形時，由于液滴上下兩端曲率的增加，表面張力也隨之增加，與磁界面力趨于平衡，鐵磁液滴的形態也達到了平衡態。隨著磁場強度的增大，鐵磁流體液滴將會在垂直方向上被進一步被拉長。圖11（b）將求得的平衡態的長短軸比 b/a 的數值結果與實驗結果[3]和 Hu 等\"的數值結果進行了對比，結果表明本文提出的方法能給出與實驗和參考文獻定量一致的結果。

2.4斯托克斯流動中液滴在靜態磁場作用下的變形

外加磁場能夠顯著影響鐵磁流體液滴的變形，從而能夠控制簡單剪切流動中液滴界面變形的控制，算例設置如圖12所示。半徑為 R 的液滴放置在 16R×8R 的計算域的中心，流場區域采用 200×100 的網格分辨率，計算域上下固壁以速度值 U 沿相反方向運動，左右邊界使用周期性邊界條件。磁場方向與流場方向之間的夾角為a ，引入磁鍵數 反映磁場強度。選取參數為 R e=1 ， C a=0.1 ， λ=1 ， ζ=2 ，

2.4.1 時磁流體液滴剪切變形

當施加 方向均勻磁場時，圖13展示了不同磁鍵數 下鐵磁流體液滴界面形態隨時間的演變。如圖13（a）所示，當無磁場分布時（ ），展示了鐵磁流體液滴在剪切流動下的變形過程。在慣性力、黏性力和表面張力的共同作用下，液滴向流場方向的拉伸，在 t=5 s時液滴達到平衡態形狀。需要指出的是，磁場強度為零時，液滴伸長變形方向與外加磁場方向趨近一致，將不考慮磁場作用的液滴變形稱為水動力學變形。當施加一個較小的磁場強度（ B_om=2.010 6 后，如圖13（b）所示，由于磁界面力作用，鐵磁流體液滴變形加劇，但此時液滴還沒有破裂，在 t=10 s后達到平衡態形狀。繼續增大磁場強度（ ），從圖13（c）中可以看出，顯著增大的磁場作用導致液滴快速拉長并逐漸破碎成2個子液滴。顯然，當外加磁場方向與液滴水動力學變形方向一致時，磁場作用會加劇液滴變形，持續增加磁場強度會導致液滴分裂。

2.4.2 時磁流體液滴剪切變形

當外加磁場方向轉變為 ，即與液滴水動力學變形方向垂直。如圖14所示，磁場強度為零時，液滴依然只發生水動力學變形，平衡態形狀與圖13（a）相同。施加與前述算例強度相同的磁場，即 和 ，磁界面力的作用直接改變了液滴的變形方向，使之朝著 方向拉長變薄，如圖14所示，在不同磁場強度下，液滴的拉伸變形程度會隨著磁場強度的增大而增大。

為了更全面地理解外加磁場對液滴剪切變形的影響，圖15給出了較小磁場強度作用下的平衡態液滴界面形態。在較小強度的磁場作用下， 時，液滴的變形方向發生了變化，伸長方向向磁場作用方向發生了偏轉。進一步增加磁場強度， 時，液滴變形方向則繼續向磁場作用方向偏轉。顯然，低強度的磁場作用可以改變液滴的變形方向。

3結束語

采用基于同位網格的有限體積法來求解耦合的不可壓縮Navier-Stokes方程、Cahn-Hilliard方程和Maxwell方程組成的控制方程組，從而發展了含鐵磁流體兩相流的相場方法，通過對Cahn-Hillard方程的分解和高階對流格式的應用，提高了方法的界面捕捉精度和魯棒性。Zalesak's圓盤問題的數值模擬結果表明，本文方法界面捕捉精度遠高于文獻報道的二階數值方法，與Xiao等提出的基于4階譜元離散的高精度相場方法精度相當。對液滴剪切變形問題的研究表明，小雷諾數 （R e） 條件下本文方法得出的結果與Cox理論預測更為吻合。當進一步提高 R e ，本文方法成功捕捉到了液滴分離的過程，并且觀察到了4個衛星液滴，這反映了本文方法在捕捉細小界面形變行為方面具有優秀的能力。

通過精確預測均勻磁場作用下鐵磁流體液滴的變形，驗證了本文方法在描述外加均勻磁場作用的能力。經過充分的驗證后，將本文方法應用于研究均勻磁場作用下鐵磁流體液滴的剪切破碎現象。當 C a 數較小（ C a=0.1 ）且忽略外加磁場作用時，液滴不會發生破碎。當施加 方向不同強度的均勻磁場時，可以觀察到在較小磁鍵數（ ）時，鐵磁流體液滴變形程度增強，液滴進一步被拉長，但還沒有破裂；當繼續增加磁場強度（ ），可以觀察到母液滴最終分裂成2個子液滴。改變磁場方向（ ），在較小磁場強度條件下液滴的變形方向會朝磁場作用方向偏轉，而在較大磁場強度下，液滴則會完全沿著磁場方向變形，并且隨著磁場強度的持續增加，液滴的變形將會越來越劇烈。

參考文獻

[1]ZhangSTNiuXDLiQPetal.Anumericalinvestigationthefatinofferofluiddroplets[]hysicsf Fluids， 2023，35（1）： 012102.

[2]Khan A，NiuXD，LiYetal.otiondeformationandcoalescenceoffeofluiddropletssubjectedtoaunifomagnetic field[J]. International Journal for Numerical Methods in Fluids， 2020， 92（11）： 1584-1603.

[3]WuZGTrollJJong HHetal.Aswarmofslipprymicropropelers penetratesthevitreusodyoftheeye[]cience Advances， 2018， 4（11）： eaat4388.

[4]FaXJSunMM，SunLNetal.Ferofluiddropletsasliquidmicrorobotswithmultipledeformabilities[]Advanced Functional Materials， 2020， 30（24）： 2000138.

[5]Medina-Sanchez M， MagdanzV， Guix M，et al.Swimming microrobots：soft，reconfigurable， andsmart[]Advanced Functional Materials， 2018， 28（25）： 1707228.

[6]SunMMHaoB，Yang SH，etal.Exploitigferofluidicwetting forminiature softmachines[]NatureCommunications， 2022，13：7919.

[7GaoDH，MorleyNBDhirV.Understanding magnetifieldgradienteffect fromaliquidmetaldropletmovement[]ual of Fluids Engineering，2004，126（1）：120-124.

[8]KorlieMSMuhejeeA，NitaBGetal.Modelingbubblesanddroplets inmagneticfluids[]oualofPhysicsCdene Matter： an Institute of Physics Journal， 2008， 20（20）： 204143. L9Cunnaucasnriqueravanveaaygoarartare-nuceacooeruemusoeoge break-up in shear flows[J]. Physics of Fluids， 2018， 30（12）： 122110.

[10]NochettoRH，SalgadoAJ，Tomas I.Adiffuse interface modelfortwo-phase ferrofluidflows[].ComputerMethods in Applied Mechanics and Engineering，2016， 309： 497-531.

[11]HuY，LiDC，NiuXD.hase-field-basedlatticeBoltzannmodelfor multiphase ferrofluidflows[]hysicalReviewE， 2018，98（3）： 033301.

[12]KanA，NiuXD，LiQZetal.Dynamic studyofferrodropletandbubblesmerging inferofluidbyasimplifiedmultipase lattice Boltzmann method[J]. Journal of Magnetism and Magnetic Materials， 2020， 495： 165869.

[13]Tag TQiaoZH.Effcient numericalmethods forphase-fieldequations[]ScientiaSinicaMathematica，22（6）：775.

[14]uPTFengLiuCetalAifusetefaceethdfmulatingwhaseflowsfplexfluds[]ufFluid Mechanics，2004， 515： 293-317.

[15]Dong S，Shen J.Atime-stepping scheme involvingconstantcoeficientmatricesforphase-fieldsimulationsof two-phase incompressble flows withlargedensityratios[J].JournalofComputational Physics，2012，231（17）：5788-5804.

[16]Jiag GS，ShuCWEicientimplementationofweightedENOschemes[].JouralofComputationalPhysics，9966（1）： 202-228.

[17]uang ZY，LinG，ArdekaniAM.Consistentessentiallyconservativeandbalanced-forcephase-fieldmethodtomodel incompressible two-phase flows[J]. Journal of Computational Physics， 2020， 406： 109192.

[18]Guermond JL，MinevP，Shen J.Anoverviewof projectionmethodsforincompressbleflows[].Computer Methods in Applied Mechanics and Engineering， 2006， 195（44/45/46/47）： 6011-6045.

[19]RheCM，ChowWL.Numericalstudyoftheturbulentflow past anairfoilwithtrailingedge separation[]AIAA Joual， 1983，21（11）： 1525-1532.

[20]uang ZY，Lin GArdekaniAM.Amixedupwind/centralWENOschemeforincompresible twohase flows[]Jouof Computational Physics， 2019， 387（C）： 455-480.

[21]Kim J.A generalizedcontinuous surface tension forceformulationforphase-fieldmodels for multi-component immiscible fluid flows[J].Computer Methods in AppliedMechanics and Engineering，2009，198（37/38/39/40）：3105-3112.

[22]LiuHH，ValocchiAJZhang YHetal.Laiceoltzmannphase-fieldmodelingofthermcapillryflows inacfied microchannel[J]. Journal of Computational Physics， 2014， 256（C）： 334-356.

[23]Wang HL，ChaiZHShiBC，etal.Comparative studyofthelaticeBoltzmannmodels forAlen-CahnandCahniliard equations[J]. Physical Review E， 2016， 94（3）： 033304.

[24]Liang H，ShiBC， Guo ZLetal.Phase-field-basedmultiple-relaxation-timelatice Boltzmannmodelforincompressble multiphase flows[J].Physical Review E，Statistical， Nonlinear，and SoftMaterPhysics，2014，89（5）：053320.

[25]XiaoYZeng ZZhangLQetalAhighlyaccuratebund-preservinghasefieldmethodforincompressblewohaseflows [J]. Physics of Fluids， 2022， 34（9）： 092103.

[26]TaylorGITheviscosityofafluidcontainingsmalldropsofanotherfluid[]ProceedingsoftheRoyalSocietyof Londn Series A， Containing Papers of a Mathematical and Physical Character， 1932， 138（834）： 41-48.

[27]GuidoSillner-imensalsapeofadndersipleseaflow[]Jualofhology91.

[28]KennedyMR，PozrikidisC，Skalak R.Motionanddeformationofiquiddrops，andtheheologyofdiluteemulsiosin simple shear flow[J]. Computers amp; Fluids， 1994， 23（2）： 251-278.

[29]CoxRG.Thedeformationofadropina general time-dependentfluid flow[].JournalofFluid Mechanics，1969，37（3）： 601-623.

[30]Hassan MR，WangC.Magneticfieldinduced ferofluid dropletbreakupina simple shearflowatalow Reynoldsnumber[]. Physics of Fluids， 2019， 31（12）： 127104.

[31]HuY，LiDC，JinLCetalHybridAllen-Cahn-basedlattceBoltzmannmodelforincompressbletwo-phase flowsthe reduction of numerical dispersion[J]. Physical Review E， 2019， 99（2）： 023302.

[32]LiQZ，LuZL，Chen Zetal.Aneffcientsimplified phase-fieldlaticeBoltzmannmethodforsuper-large-density-ratio multiphase flow[J]. International Journal of Multiphase Flow， 2023， 160： 104368.

[33]FlamentCLacis SBacrJCetalMeasurementsofferoflidsufaceensioinofinedgeometr]hysicalReviewE， 1996， 53（5）： 4801-4806.

（編輯鄭潔）