基于<bold>WENO</bold>格式有限体积法的铁磁流体两相流相场方法

张少松，张良奇，陈黎明，王小双，肖姚，曾忠; ZHANG Shaosong; ZHANG Liangqi; CHENG Liming; WANG Xiaoshuang; XIAO Yao; ZENG Zhong

网刊加载中。。。

使用Chrome浏览器效果最佳，继续浏览，你可能不会看到最佳的展示效果，

确定继续浏览么?

复制成功，请在其他浏览器进行阅读

基于WENO格式有限体积法的铁磁流体两相流相场方法 PDF

- ORCID：
张少松 ¹
✉
- ORCID：
张良奇 ¹
✉
- ORCID：
陈黎明 ¹
- ORCID：
王小双 ²
- ORCID：
肖姚 ¹
- ORCID：
曾忠 ¹

1. 重庆大学航空航天学院，重庆 400044； 2. 沪渝人工智能研究院，重庆 401332

中图分类号： O359+.1

最近更新：2025-04-25

DOI：10.11835/j.issn.1000-582X.2024.271

摘要

采用基于WENO格式的有限体积法，发展了包括铁磁流体的两相流相场方法。采用不可压缩Navier-Stokes方程描述流体流动，利用Cahn-Hilliard方程捕捉两相流界面运动，并采用Maxwell方程描述磁场分布。同时，在流体流动控制方程中加入开尔文力和表面张力以实现磁场对界面动力学行为的描述。将四阶Cahn-Hilliard方程拆分为2个Helmholtz类型方程，从而克服四阶非线性项的离散和高精度计算带来的难题。采用五阶WENO格式对控制方程的对流项进行统一离散处理，从而提高数值计算的精度性，同时避免产生数值振荡。Zalesak’s圆盘问题数值模拟的结果表明，本文方法的相界面捕捉精度高于参考文献中所报道的离散方法，与高精度相场方法的精度相当；对剪切流中的液滴变形问题的数值模拟揭示本文数值方法可以更真实地捕捉到更多个卫星液滴。此外，针对磁场影响下，较低毛细数条件的铁磁流体液滴剪切变形研究显示，当外加磁场方向与液滴水动力学变形方向一致时，磁场的作用会放大液滴变形，进一步增加磁场强度会诱发液滴分裂；而当外加磁场方向与液滴水动力学变形方向垂直时，较低强度的磁场作用能够改变液滴变形方向，而较高强度的磁场则会使液滴直接呈现出沿磁场方向变形。

关键词

相场法; 有限体积法; 铁磁流体液滴; 磁场控制

铁磁流体凭借其优越的变形能力和磁场可控性而受到广泛应用^[

1]。最具代表性的应用是铁磁流体液滴可用作反应平台或载体，在芯片实验室中实现高效可控的样品输运。此外，铁磁流体液滴在外加磁场的精确控制下能分裂成多个子液滴，再按需重新融合，被视为能实现复杂功能的微型磁性软机器人，在医疗保健、生物工程和微操作等领域展现了巨大应用潜力，吸引了大量研究者的兴趣^{[参考文献 2⁃6}2⁃6]。

为了深入研究铁磁流体的界面动力学行为，已经发展出多种多相铁磁流体流动数值方法。Gao等^[

7]采用流体体积（VOF）方法模拟了微重力真空条件下金属液滴在磁场区域内的运动情况。Korlie等^{[参考文献 8

百度学术}8]提出了一种VOF方法来模拟铁磁流体中气泡和液滴的运动。Cunha Lucas等^{[参考文献 9

百度学术}9]利用Level Set方法对在较高雷诺数（Re=1）下，磁场诱导的剪切流中铁磁液滴破碎的现象进行数值模拟研究。研究结果表明，当磁场作用方向与流动方向平行时，会延迟破裂过程，减小卫星液滴的尺寸；相反，如果在垂直于流动方向上施加磁场，则可以通过调整磁场强度来诱导、延迟甚至阻止液滴的破裂。近期有许多学者开始使用相场法来描述多相铁磁流体流动行为。与前述的2种界面描述方法相比，相场法避免了额外的界面几何重构和重新初始化过程，而且表面张力计算简单，并且能够确保质量和能量的守恒。Nochetto等^{[参考文献 10

百度学术}10]开发出含铁磁流体两相流的相场模型，能捕获铁磁流体的基本流动现象，如Rosensweig不稳定性和铁磁流体刺猬。Hu等^{[参考文献 11

百度学术}11]将相场方法应用于铁磁流体流动，采用守恒的Allen-Cahn方程捕捉界面。他们开发了一种多松弛格子玻尔兹曼方法（MRT-LBM）来求解控制方程。Khan等^{[参考文献 12

百度学术}12]用新近发展的简化多相格子玻尔兹曼方法（SMLBM）捕捉铁磁流体的流体动力学行为，通过数值模拟研究，探讨了在均匀磁场作用下铁磁流体中铁磁液滴变形和磁流体气泡融合的动力学过程。

铁磁流体液滴在微流控技术和生物医药领域的前沿应用需要精确控制其界面动力学行为，而对铁磁液滴运动和变形行为进行数值模拟研究能为理解控制机理和设计控制方式提供必要的理论支持。笔者拟选择铁磁流体相场模型为理论模型，应用基于五阶WENO格式的有限体积方法统一离散求解控制方程组，以准确描述铁磁流体液滴的界面动力学行为。Tang等^[

13]指出，基于Cahn-Hilliard方程相场模型的数值求解存在很大难度，如处理高阶非线性扩散项、保持长时间数值计算的精度，以及能量稳定性等，都是需要克服的挑战。Yue等^{[参考文献 14

百度学术}14]明确指出相变量的数值扰动会导致液滴体积的非物理收缩，并破坏质量守恒。这里，借鉴了Dong等^{[参考文献 15

百度学术}15]的工作思路，将四阶Cahn-Hilliard方程拆分为2个Helmholtz形式方程，这样可以克服由四阶非线性项带来的困难。为了提升鲁棒性，采用了五阶WENO格式来统一处理控制方程的对流项，以提高数值计算的精度，同时避免数值振荡^{[参考文献 16

百度学术}16]。此外，参考Huang等^{[参考文献 17

百度学术}17]的研究，引入了一致、守恒的质量通量来描述相变量扩散对质量和动量输运的影响，保证相界面处动量输运的守恒性。采用增量压力修正方案，实现不可压缩流的流动速度和压力的解耦，从而在每一个时间步，只需要求解一系列解耦的速度和压力离散方程^{[参考文献 18

百度学术}18]。本文提出的方法采用了统一的控制体构造方式，所有变量都存储在控制体中心，也就是采用同位网格。采用了Rhie-Chow插值修正面心速度，以抑制压力振荡^{[参考文献 19

百度学术}19]。在采用Zalesak’s圆盘问题验证所提出方法在捕捉界面精度上的优越性后，将其应用于研究在磁场作用下铁磁流体液滴剪切变形的动力学。

1 理论模型与数值方法

1.1 两相铁磁流体相场模型

研究对象是包含铁磁流体的两相流，其中一相为铁磁流体，另一相为非磁性流体。两相铁磁流体流动的相场模型归纳如下。

1.1.1 相场方程

采用Cahn-Hilliard相场模型^[

20]来描述不可压缩两相流的非混溶界面。

\frac{\partial ϕ}{\partial t} + \nabla \cdot (u ϕ) = \nabla \cdot (M \nabla ξ)

，

（1）

ξ = K (F' (ϕ) - \nabla^{2} ϕ)

，

（2）

K = \frac{3}{2 \sqrt[]{2}} σ ε

，

（3）

F (ϕ) = \frac{1}{4 ε^{2}} {(1 - ϕ)}^{2} {(1 + ϕ)}^{2}

，

（4）

式中： $ϕ = 1$ 和 $- 1$ 分别为流体1和流体2；u为速度；M为迁移率；ξ为化学势；K为混合能密度；ε为相场模型中的界面厚度；σ为2种流体之间的表面张力系数；F( $ϕ$ )为体自由能密度，通常选用一个双阱势函数来表示。

1.1.2 流动方程

采用不可压缩Navier-Stokes方程描述两相流。同时，参考了Huang等^[

17]的研究工作，引入一致、守恒的质量通量m_cc来描述相变量扩散对质量、动量输运的影响，以确保相界面处动量输运的守恒性。

\frac{\partial ρ u}{\partial t} + \nabla \cdot (m_{c c} \otimes u) = - \nabla p + \nabla \cdot [μ (\nabla u + {(\nabla u)}^{T})] + F_{t o t a l}

，

（5）

\nabla \cdot u = 0

，

（6）

m_{c c} = ρ u - \frac{ρ_{1} - ρ_{2}}{2} M \nabla ξ

，

（7）

式中：ρ是密度；μ是动力黏度；p是压力；F_total是外力项（包括表面张力F_σ和开尔文力F_mg）。

应用连续表面力模型^[

21⁃22]计算表面张力。当流场的表面张力系数为常数时，表面张力则满足

F_{σ} = \frac{3 \sqrt[]{2}}{4} \frac{σ}{ε} ξ \nabla ϕ

。

（8）

两相体系的密度ρ和动力黏度μ与序参量 $ϕ$ 之间存在线性关系，定义为

ρ = \frac{ρ_{1} + ρ_{2}}{2} + \frac{ρ_{1} - ρ_{2}}{2} ϕ, μ = \frac{μ_{1} + μ_{2}}{2} + \frac{μ_{1} - μ_{2}}{2} ϕ

，

（9）

式中：ρ₁和ρ₂分别为2种流体的密度；μ₁和μ₂分别为2种流体的动力黏度。

1.1.3 磁场方程

由于关注重点在铁磁流体的界面动力学，故忽略了铁磁流体的内旋效应，而采用非导电铁磁流体的Maxwell方程^[

11]来描述磁场。

\nabla \cdot B = 0

，

（10）

\nabla \cdot Η = 0

，

（11）

式中：B和H分别为磁感应强度和磁场强度。B可以表示为

B = η_{0} (H + M) = η_{0} (1 + χ) H = η H

，

（12）

式中：M为磁化强度；η是磁导率；χ是磁化率；η₀为真空磁导率，这里取值为4π×10^-7 N/A²。根据方程（11）所示磁场H的无旋条件，引入一个磁标量势 $ψ$ ，定义为

H = - \nabla ψ

。

（13）

将方程（13）代入方程（12），然后代入方程（10），得到磁势方程

\nabla \cdot (η \nabla ψ) = 0

。

（14）

在计算过程中，通过更新磁势方程（14）来求解磁标量势 $ψ$ 。根据方程（17），磁标量势 $ψ$ 与相界面 $ϕ$ 之间通过磁导率η耦合。当相界面发生变化时，对应的磁势也会随之改变，从而更新磁场H。

外加磁场会对铁磁流体相产生附加磁应力

τ_{m} = - \frac{η_{0}}{2} | H |^{2} I + H B

，

（15）

式中：τ_m表示附加的磁应力张量，磁应力张量的作用等价于在磁流体相中施加开尔文力F_mg。

F_{m g} = \nabla \cdot τ_{m} = - \frac{η_{0}}{2} \nabla (| H |^{2}) + (B \cdot \nabla) H = - \frac{η_{0} χ}{2} \nabla | H |^{2}

，

（16）

两相体系的磁导率 $η$ 和磁化率 $χ$ 与序参量 $ϕ$ 之间存在线性关系

η = \frac{η_{1} + η_{2}}{2} + \frac{η_{1} - η_{2}}{2} ϕ, χ = \frac{χ_{1} + χ_{2}}{2} + \frac{χ_{1} - χ_{2}}{2} ϕ

，

（17）

式中：η₁和η₂分别为2种流体的磁导率；χ₁和χ₂分别为2种流体的磁化率。

垂直磁场情况下，上下壁面边界条件为 $\frac{\partial ψ}{\partial y} = H$ ，左右壁面边界条件为 $\frac{\partial ψ}{\partial x} = 0$ ， $H$ 为初始磁场强度。

1.2 控制方程的离散

采用同位网格进行离散，即采用统一的控制体构造方式，而且所有变量都存储在控制中心，对流项采用五阶WENO显式处理。

1.1.1 时间离散

瞬态项采用二阶向后差分进行离散，

\frac{\partial f}{\partial t} = \frac{3 f^{n + 1} - 4 f^{n} + f^{n - 1}}{2 Δ t}

，

（18）

式中：f为标量函数的节点值； $Δ t$ 为时间步长；用上标n表示当前时间步；n+1表示下一时间步； $n - 1$ 表示上一时间步。

为了方便后续离散的描述，瞬态项离散可以写为

\frac{\partial f^{n + 1}}{\partial t} = \frac{γ_{0} f^{n + 1} - \overset{⌢}{f}}{Δ t}

，

（19）

其中

γ_{0} = \{\begin{array}{l} 1, & i f J = 1 ， \\ 1.5, & i f J = 2 ， \end{array}

（20）

\overset{⌢}{f} = \{\begin{array}{l} f^{n}, & i f J = 1, \\ 2 f^{n} - \frac{1}{2} f^{n - 1}, & i f J = 2 。 \end{array}

（21）

除此之外，定义 $f^{*, n + 1}$ 代表时间上J阶显式逼近 $f^{n + 1}$ ，

f^{*, n + 1} = \{\begin{array}{l} f^{n}, & i f J = 1 ， \\ 2 f^{n} - f^{n - 1}, & i f J = 2 。 \end{array}

（22）

1.1.2 空间算子离散

基于中心差分格式，对控制方程的梯度、散度、拉普拉斯算子进行离散，梯度算子离散为

\nabla f (i, j) = [\begin{array}{l} \partial f / \partial x \\ \partial f / \partial y \end{array}] = [\begin{array}{l} \frac{f (i + 1, j) - f (i - 1, j)}{2 Δ x} \\ \frac{f (i, j + 1) - f (i, j - 1)}{2 Δ y} \end{array}]

，

（23）

式中： $Δ x$ 和 $Δ y$ 分别表示x轴和y轴上的网格步长。

散度算子离散为

\nabla \cdot F (i, j) = \frac{\partial F_{x}}{\partial x} + \frac{\partial F_{y}}{\partial y} = \frac{F_{x} (i + 1, j) - F_{x} (i - 1, j)}{2 Δ x} + \frac{F_{y} (i, j + 1) - F_{y} (i, j - 1)}{2 Δ y}

，

（24）

式中： $F (i, j)$ 为矢量函数的节点值；F_x为矢量函数x轴上的节点值；F_y为矢量函数y轴上的节点值。

拉普拉斯算子离散为

\begin{array}{l} \nabla^{2} f (i, j) = \nabla \cdot \nabla f (i, j) = \\ \frac{\partial_{x} f (i + \frac{1}{2}, j) - \partial_{x} f (i - \frac{1}{2}, j)}{Δ x} + \frac{\partial_{y} f (i, j + \frac{1}{2}) - \partial_{y} f (i, j - \frac{1}{2})}{Δ y} = \\ \frac{\frac{f (i + 1, j) - f (i, j)}{Δ x} - \frac{f (i, j) - f (i - 1, j)}{Δ x}}{Δ x} + \frac{\frac{f (i, j + 1) - f (i, j)}{Δ y} - \frac{f (i, j) - f (i, j - 1)}{Δ y}}{Δ y} = \\ - (\frac{2}{Δ x^{2}} + \frac{2}{Δ y^{2}}) f (i, j) + \frac{1}{Δ x^{2}} f (i - 1, j) + \frac{1}{Δ x^{2}} f (i + 1, j) + \frac{1}{Δ y^{2}} f (i, j - 1) + \frac{1}{Δ y^{2}} f (i, j + 1) 。 \end{array}

（25）

1.1.3 Cahn-Hilliard方程的离散

将相场方程分解为2个二阶Helmholtz方程离散求解^[

15]

\{\begin{array}{l} \nabla^{2} φ - (α + \frac{S}{ε^{2}}) φ = Q ， \\ \nabla^{2} ϕ^{n + 1} + α ϕ^{n + 1} = φ 。 \end{array}

（26）

其中

S \geq ε^{2} \sqrt[]{\frac{4 γ_{0}}{K M Δ t}}

，

（27）

α = - \frac{S}{2 ε^{2}} (1 + \sqrt[]{1 - \frac{4 γ_{0}}{K M Δ t} \frac{ε^{4}}{S^{2}}})

，

（28）

Q = \frac{1}{K M} [\frac{\hat{ϕ}}{Δ t} - \nabla \cdot (u^{*, n + 1} {\tilde{ϕ}}^{*, n + 1}) + S_{ϕ}^{n + 1}] + \nabla^{2} F' (ϕ^{*, n + 1}) - \frac{S}{ε^{2}} \nabla^{2} ϕ^{*, n + 1}

，

（29）

式中： $γ_{0}$ 值为1.5，当时间格式为二阶向后差分时， $\hat{(\cdot)} = 2 {(\cdot)}^{n} - 0.5 {(\cdot)}^{n - 1}$ ， ${(\cdot)}^{*, n + 1} = 2 {(\cdot)}^{n} - {(\cdot)}^{n - 1}$ ，对流项中的变量插值 $\tilde{(\cdot)}$ 采用五阶WENO格式处理，其他空间算子均用中心差分格式。 $ϕ^{n + 1}$ 和 $φ$ 的边界条件为

n \cdot \nabla ϕ = 0, n \cdot \nabla φ = 0

。

（30）

1.1.4 Navier-Stokes方程的离散

采用增量压力修正投影法对Navier-Stokes方程进行解耦求解，将动量方程（5）拆分为方程（31）和方程（32）。

\frac{\partial (ρ u)}{\partial t} |_{n}^{*} + \nabla \cdot (m_{c c} \otimes u) = \nabla \cdot [μ (\nabla u + \nabla u^{T})] - \nabla p^{n} + F_{t o t a l}

，

（31）

\frac{\partial (ρ u)}{\partial t} |_{*}^{n + 1} = \nabla (- p^{n + 1} + p^{n})

，

（32）

式中：u^*为中间速度。

1）预测步。将方程（31）的部分黏性项移到方程左侧进行隐式处理，对流项采用五阶WENO显式处理，如下所示。

\begin{array}{l} \frac{\frac{3}{2} ρ^{n + 1} u^{*}}{d t} - \nabla \cdot (μ^{n + 1} \nabla u^{*}) = \frac{2 ρ^{n} u^{n} - \frac{1}{2} ρ^{n - 1} u^{n - 1}}{d t} - \nabla \cdot {(m_{c c} \otimes \tilde{u})}^{*, n + 1} + \\ \nabla \cdot (μ^{n + 1} \nabla {(u^{*, n + 1})}^{T}) - \nabla p^{*, n + 1} + F_{t o t a l}^{^{n + 1}} 。 \end{array}

（33）

2）压力修正步。在方程（32）两端同时取散度，结合连续性方程（6）推导出增量压力泊松方程

\nabla \cdot (\frac{1}{ρ^{n + 1}} \nabla p^{'}) = \frac{3}{2 d t} \nabla \cdot {\hat{u}}^{*}

，

（34）

p^{'} = p^{n + 1} - p^{n}

。

（35）

3）速度插值。采用同位网格进行离散，即采用统一的控制体构造方式，而且所有变量都存储在控制中心，离散方程（34）时引入Rhie-Chow插值，以避免压力震荡。将式（33）在如图1所示的控制体积P中进行离散，得到的各项系数为

\begin{array}{l} a_{P} = \frac{μ_{W} + 2 μ_{P} + μ_{E}}{2} \frac{Δ y}{Δ x} + \frac{μ_{N} + 2 μ_{P} + μ_{S}}{2} \frac{Δ x}{Δ y} ， \\ a_{E} = \frac{μ_{E} + μ_{P}}{2} \frac{Δ y}{Δ x}, a_{W} = \frac{μ_{W} + μ_{P}}{2} \frac{Δ y}{Δ x} ， \\ a_{N} = \frac{μ_{N} + μ_{P}}{2} \frac{Δ x}{Δ y}, a_{S} = \frac{μ_{S} + μ_{P}}{2} \frac{Δ x}{Δ y} 。 \end{array}

（36）

图1 等距同位网格

Fig. 1 Equidistant parity grid

f为网格单元的面心

采用Rhie-Chow插值修正面心速度 $u_{f}^{*}$ 以实现对方程（34）右端速度散度项的处理

{\hat{u}}_{f_{x}}^{*} = \frac{u_{P}^{*} + u_{E}^{*}}{2} + R_{x}

，

（37）

{\hat{v}}_{_{f_{y}}}^{*} = \frac{v_{P}^{*} + v_{N}^{*}}{2} + R_{y}

，

（38）

式中：下标 $f_{x}$ 表示节点P和E的中间界面； $f_{y}$ 表示P和N节点的中间界面。在相场两相流中，修正算子R计算如下，

R_{x} = - \frac{{\hat{d}}_{f_{x}}}{4 Δ x} [p_{W} - 3 p_{P} + 3 p_{E} - p_{E E}] + \frac{3 \sqrt[]{2}}{4} \frac{σ}{ε} [ξ_{f_{x}} \frac{\partial ϕ}{\partial x} |_{f_{x}} - \frac{ρ_{f_{x}}}{2} (\frac{ξ_{P}}{ρ_{P}} \frac{\partial ϕ}{\partial x} |_{P} + \frac{ξ_{E}}{ρ_{E}} \frac{\partial ϕ}{\partial x} |_{E})]

，

（39）

R_{y} = - \frac{{\hat{d}}_{f_{y}}}{4 Δ y} [p_{S} - 3 p_{P} + 3 p_{N} - p_{N N}] + \frac{3 \sqrt[]{2}}{4} \frac{σ}{ε} [ξ_{f_{y}} \frac{\partial ϕ}{\partial y} |_{f_{y}} - \frac{ρ_{f_{y}}}{2} (\frac{ξ_{P}}{ρ_{P}} \frac{\partial ϕ}{\partial y} |_{P} + \frac{ξ_{N}}{ρ_{N}} \frac{\partial ϕ}{\partial y} |_{N})]

，

（40）

{\hat{d}}_{f} = \frac{\frac{3}{2} \frac{ρ_{f}}{d t} \cdot \frac{Δ V}{a_{P}}}{1 + \frac{3}{2} \frac{ρ_{f}}{d t} \cdot \frac{Δ V}{a_{P}}}

，

（41）

式中： $f$ 表示面心； $Δ V$ 为单元的体积。

4）最终步。由前两步求得的 ${\hat{u}}^{*}$ 和 $p^{'}$ 求解出 $u^{n + 1}$ 。

u^{n + 1} = {\hat{u}}^{*} - \frac{2 d t}{3 ρ^{n + 1}} \nabla p^{'}

。

（42）

2 数值验证

2.1 界面捕捉精度的数值验证

利用经典的Zalesak’s圆盘旋转案例来验证本文方法在界面捕捉精度方面的效果。理论上，经过一个或者多个周期的旋转后，圆盘将回到初始位置，因此，通过比较经过周期旋转后的界面形态与初始界面形态即可反映数值计算界面捕捉精度。算例设置如图2所示。一个带凹槽的圆盘放在大小为L₀×L₀的周期区域的中心，取L₀=200，圆盘的半径和凹槽的宽度分别为80和16，扩散界面区域的厚度设置为ε=2.0。圆盘在如下速度场的驱动下做周期性旋转。

u = - U_{0} π (\frac{y}{L_{0}} - 0.5), v = U_{0} π (\frac{x}{L_{0}} - 0.5)

。

（43）

图2 Zalesak’s 圆盘的初始形状

Fig. 2 The initial shape of Zalesak’s disk

实际计算中，取特征速度U₀=0.02和0.04。时间步长和网格步长为∆t=∆x=1，四周采用纽曼边界条件。定义无量纲参数Peclet´为Pe=U₀ε/M以反映算例设置。

图3展示了在Pe=800时的圆盘周期性旋转的过程。从图3可以看出在经历一个完整的周期后，圆盘的形状与初始形状几乎完全重合。此外，为进一步量化当前方法的界面捕捉精度，定义如下误差以定量表征数值结果的精度。

图3 Pe=800时不同时刻的界面形状

Fig. 3 When Pe=800, the interface shapes at different time points

E_{r} = \frac{\sum_{i} |ϕ_{_{i}}^{f i n a l} - ϕ_{_{i}}^{0}|}{\sum_{i} |ϕ_{_{i}}^{0}|}

，

（44）

式中： $ϕ_{_{i}}^{f i n a l}$ 为数值计算得到的一个周期后相变量分布； $ϕ_{_{i}}^{0}$ 为初始相变量分布。表1对比分析了Pe=400和800下界面计算误差。在相同的网格分辨率下，当前方法数值计算的误差远小于文献报道的结果^[

23⁃24]：对于2个M=0.000 1的算例，当前计算结果的误差约为Wang等^{[参考文献 23

百度学术}23]报道误差的1/20。本文方法精度与Xiao等^{[参考文献 25

百度学术}25]提出的基于4阶谱元离散的高精度相场方法相当，这充分反映了当前方法的良好界面捕捉精度。

表1 U₀=0.02和0.04时，Zalesak圆盘问题本文方法与文献[23⁃25]的相对误差对比

Table 1 When U₀=0.02 and 0.04, comparison of the Zalesak disk problem’s relative errors of method presented in this paper with those of references [23-25]

M	U₀=0.02				U₀=0.04
M	本文方法	文献[23]（AC）	文献[24]（CH）	文献[25]	本文方法	文献[23]（AC）	文献[24]（CH）	文献[25]
0.01	0.016 507 0	0.026 958	0.042 114	0.025 6	0.010 937	0.026 641	0.048 562	0.012 6
0.001	0.004 900 9	0.041 451	0.055 458	0.007 8	0.003 567 4	0.040 926	0.069 678	0.003 8
0.000 1	0.001 732 8	0.043 452	0.095 332	0.002 3	0.002 712 8	0.049 172	—	0.001 8

图4展示了Pe=800时，不同网格分辨率N×N的误差变化曲线，横坐标为网格分辨率N（N=50、100、200、256），纵坐标为误差E_r。结果表明，误差随网格分辨率呈线性收敛，收敛阶数约为2.5。

图4 误差随网格收敛的变化曲线图

Fig. 4 Error variation with grid convergence plot

2.2 液滴在剪切流动下的变形

单个液滴在剪切流动中的变形是检验两相流数值方法界面捕捉精度的另一个经典算例，该算例设置如图5所示。一个半径为R的液滴最初放置在8R×8R的计算域中心，网格分辨率为200×200。计算域的上下固壁分别以大小为U的速度向相反方向移动，左右边界为周期性边界条件。为了量化液滴变形，定义液滴的变形参数 $D = (L - B) / (L + B)$ ，L和B分别为变形液滴在平衡状态下的长轴和短轴的长度。

图5 剪切流动中液滴变形的初始形状

Fig. 5 The initial deformation of droplets in shear flow

当前算例界面变形的无量纲控制参数有

C a = μ \dot{γ} R / σ

，

R e = ρ \dot{γ} R^{2} / μ

，

（45）

式中：Ca为毛细数；Re为雷诺数； $\dot{γ}$ 为剪切率；λ为黏度比。

λ = \frac{μ_{h}}{μ_{l}}, \dot{γ} = U / (L_{y})

。

（46）

文献[

26]指出，当Ca<<1和Re<<1时，液滴剪切变形满足如下理论关系

D_{t a y l o r} = \frac{L - B}{L + B} = \frac{19 λ + 16}{16 λ + 16} C a

。

（47）

上式变形参数预测依赖于流动无界且雷诺数趋近于零的假设。然而，在数值研究中，上下2个运动壁面会对内部流动产生约束效应，定义几何约束比为2R₀/H_d以量化约束效应的影响。文献[

27⁃28]指出，当约束比2R₀/H_d<0.4时，约束效应对液滴变形的影响可以忽略不计，所以选取2R₀/H_d=0.25。Cox^{[参考文献 29

百度学术}29]对Taylor理论做了延伸拓展以描述较大Ca数条件下的液滴变形

D_{c o x} = D_{t a y l o r} {[1 + {(\frac{19 C a λ}{20})}^{2}]}^{- \frac{1}{2}}

。

（48）

图6展示了不同毛细数下（Ca=0.1、0.25、0.4、0.5）平衡态液滴界面形态，可以看出随着Ca数的不断增大，黏性力作用增强，液滴变形程度越大。图7给出了当前数值计算得到的变形参数随Ca的变化，当前方法所得数值结果与Cox变形理论和参考文献[

30]非常吻合。

图6 不同毛细数下平衡液滴形态

Fig. 6 Balanced droplet morphologies at different capillary numbers

图7 数值结果与Cox理论和Level Set方法的对比

Fig. 7 Numerical results comparison between Cox theory and Level Set method

为验证当前数值方法捕捉液滴破碎行为的能力，探讨了在更高雷诺数（Re）条件下的液滴剪切变形。取Re=1，Ca=0.5，λ=1.2，Cn=ε/L_y=0.007 5；半径为R的液滴放置在16R×8R的计算域的中心，网格分辨率为200×100。

图8展示了不同时刻液滴界面形状的演化过程，观察到液滴破裂前的界面变形结果与文献[

9]和文献[30]的报道是一致的。然而，在液滴破碎以后，本文的相场方法预测出现了4个卫星液滴，而相同条件下，文献[9]报道只出现了1个卫星液滴，文献[30]的结果出现了2个卫星液滴。Hu等^{[参考文献 31

百度学术}31]和Li等^{[参考文献 32

百度学术}32]指出，改善守恒性、减少数值耗散能够显著提高相场方法捕捉微小界面变形细节的能力，他们改良的相场方法分别在液滴剪切问题和液滴撞击液膜问题中捕捉到了更多的卫星液滴。此外，相场方法界面捕捉精度与界面厚度和空间网格分辨率相关。这里首先对当前的数值计算进行Cn数和网格分辨率的收敛性验证。如图9所示，较大Cn数对应过宽的界面厚度，使微小界面结构变化被抹平，所以在Cn=0.01时，没有观测到卫星液滴；降低Cn数，数值计算捕捉细小界面变形的能力显著改善，当Cn=0.008时出现4个卫星液滴，再进一步降低界面厚度，不再出现更多卫星液滴，所以当前结果是关于界面厚度参数收敛的。此外，如图10所示，固定Cn=0.007 5，网格分辨率增大，液滴破碎后仍观察到4个卫星液滴，所以网格分辨率为200×100时数值结果已经实现收敛。

图8 当Re=1，Ca=0.5，λ=1.2时，单液滴剪切流中液滴破碎随时间变化图

Fig. 8 Evolution of the broken of a droplet in a single droplet shear flow over time when Re=1, Ca=0.5, and λ=1.2

图9 当Re=1，Ca=0.5，λ=1.2时，单液滴剪切流中不同Cn数液滴破碎图

Fig. 9 The broken of a droplet in a single droplet shear flow with different Cn numbers when Re=1, Ca=0.5, and λ=1.2

图10 当Re=1，Ca=0.5，λ=1.2，Cn=0.007 5时，网格分辨率为400×200时液滴破碎图

Fig. 10 The broken of a droplet with mesh resolution of 400×200 when Re=1, Ca=0.5, λ=1.2, and Cn=0.007 5

2.3 静态磁场作用下液滴的变形

为验证当前方法描述磁场作用的准确性，模拟了不同磁场强度条件下的铁磁流体液滴变形。在这种情况下，将半径为1的铁磁流体液滴放置在16×16的矩形区域的中心。铁磁流体和非磁性环境流体的材料性质与Flament实验研究^[

33]类似。当施加垂直向上方向的均匀磁场时，开尔文力会在2种流体界面处由于磁场作用产生，并打破表面张力与压力梯度的平衡，铁磁流体液滴发生沿磁场方向的变形。选用的参数为黏度比λ=1，磁导率比ζ=η_h/η_l=3.2，网格分辨率为200×200。

图11（a）描述了在不同磁场强度H（H=1.2，2.4，2.9，3.7，5.5 kA/m）下，铁磁流体液滴的平衡态界面形状。当施加竖直方向的均匀磁场时，液滴会在磁界面力的作用下沿磁场方向伸长，而表面张力的作用则努力使液滴恢复到圆形。当液滴开始变形时，由于液滴上下两端曲率的增加，表面张力也随之增加，与磁界面力趋于平衡，铁磁液滴的形态也达到了平衡态。随着磁场强度的增大，铁磁流体液滴将会在垂直方向上被进一步被拉长。图11（b）将求得的平衡态的长短轴比b/a的数值结果与实验结果^[

33]和Hu等^{[参考文献 11

百度学术}11]的数值结果进行了对比，结果表明本文提出的方法能给出与实验和参考文献定量一致的结果。

图11 不同磁场强度H时磁流体液滴的平衡态形状图和长径比与实验结果^[

33]和参考文献^{[参考文献 11

百度学术}11]的比较图

Fig. 11 The equilibrium shape and comparison of the length-diameter ratio of magnetic fluid droplets with experimental results^[

33] and reference^{[参考文献 11

百度学术}11] at different magnetic field intensities H

2.4 斯托克斯流动中液滴在静态磁场作用下的变形

外加磁场能够显著影响铁磁流体液滴的变形，从而能够控制简单剪切流动中液滴界面变形的控制，算例设置如图12所示。半径为R的液滴放置在16R×8R的计算域的中心，流场区域采用200×100的网格分辨率，计算域上下固壁以速度值U沿相反方向运动，左右边界使用周期性边界条件。磁场方向与流场方向之间的夹角为α，引入磁键数B_om=Rη₀H $_{0}^{2}$ /2σ反映磁场强度。选取参数为 $R e = 1 ， C a = 0.1 ， λ = 1 ， ζ = 2 ， C n = ε / L_{y} = 0.007 5$ 。

图12 静态磁场作用下液滴在剪切流动中的初始形状

Fig. 12 Initial shape of droplet in shear flow under static magnetic field

2.4.1 α=45°时磁流体液滴剪切变形

当施加45°方向均匀磁场时，图13展示了不同磁键数B_om下铁磁流体液滴界面形态随时间的演变。如图13（a）所示，当无磁场分布时（B_om=0），展示了铁磁流体液滴在剪切流动下的变形过程。在惯性力、黏性力和表面张力的共同作用下，液滴向流场方向的拉伸，在t=5 s时液滴达到平衡态形状。需要指出的是，磁场强度为零时，液滴伸长变形方向与外加磁场方向趋近一致，将不考虑磁场作用的液滴变形称为水动力学变形。当施加一个较小的磁场强度（B_om=2.010 6）后，如图13（b）所示，由于磁界面力作用，铁磁流体液滴变形加剧，但此时液滴还没有破裂，在t=10 s后达到平衡态形状。继续增大磁场强度（B_om=12.566 4），从图13（c）中可以看出，显著增大的磁场作用导致液滴快速拉长并逐渐破碎成2个子液滴。显然，当外加磁场方向与液滴水动力学变形方向一致时，磁场作用会加剧液滴变形，持续增加磁场强度会导致液滴分裂。

图13 不同B_om数单液滴剪切流中液滴随时间变化图

Fig. 13 Droplet change over time in a single droplet shear flow with different B_om numbers

2.4.2 α=135°时磁流体液滴剪切变形

当外加磁场方向转变为 $α = 135^{o}$ ，即与液滴水动力学变形方向垂直。如图14所示，磁场强度为零时，液滴依然只发生水动力学变形，平衡态形状与图13（a）相同。施加与前述算例强度相同的磁场，即B_om=2.010 6和B_om=12.566 4，磁界面力的作用直接改变了液滴的变形方向，使之朝着135°方向拉长变薄，如图14所示，在不同磁场强度下，液滴的拉伸变形程度会随着磁场强度的增大而增大。

图14 不同B_om数单液滴剪切流中液滴平衡态图

Fig. 14 Droplet equilibrium in shear flow of single droplet with different B_om numbers

为了更全面地理解外加磁场对液滴剪切变形的影响，图15给出了较小磁场强度作用下的平衡态液滴界面形态。在较小强度的磁场作用下，B_om=0.282 7时，液滴的变形方向发生了变化，伸长方向向磁场作用方向发生了偏转。进一步增加磁场强度，B_om=0.502 7时，液滴变形方向则继续向磁场作用方向偏转。显然，低强度的磁场作用可以改变液滴的变形方向。

图15 在B_om数较小时单液滴剪切流中液滴平衡态图

Fig. 15 Droplet equilibrium in a single droplet shear flow with smaller B_om number

3 结束语

采用基于同位网格的有限体积法来求解耦合的不可压缩Navier-Stokes方程、Cahn-Hilliard方程和Maxwell方程组成的控制方程组，从而发展了含铁磁流体两相流的相场方法，通过对Cahn-Hilliard方程的分解和高阶对流格式的应用，提高了方法的界面捕捉精度和鲁棒性。Zalesak’s圆盘问题的数值模拟结果表明，本文方法界面捕捉精度远高于文献报道的二阶数值方法，与Xiao等^[

25]提出的基于4阶谱元离散的高精度相场方法精度相当。对液滴剪切变形问题的研究表明，小雷诺数（Re）条件下本文方法得出的结果与Cox理论预测更为吻合。当进一步提高Re，本文方法成功捕捉到了液滴分离的过程，并且观察到了4个卫星液滴，这反映了本文方法在捕捉细小界面形变行为方面具有优秀的能力。

通过精确预测均匀磁场作用下铁磁流体液滴的变形，验证了本文方法在描述外加均匀磁场作用的能力。经过充分的验证后，将本文方法应用于研究均匀磁场作用下铁磁流体液滴的剪切破碎现象。当Ca数较小（Ca=0.1）且忽略外加磁场作用时，液滴不会发生破碎。当施加45°方向不同强度的均匀磁场时，可以观察到在较小磁键数（B_om=2.010 6）时，铁磁流体液滴变形程度增强，液滴进一步被拉长，但还没有破裂；当继续增加磁场强度（B_om=12.566 4），可以观察到母液滴最终分裂成2个子液滴。改变磁场方向（ $α =$ 135°），在较小磁场强度条件下液滴的变形方向会朝磁场作用方向偏转，而在较大磁场强度下，液滴则会完全沿着磁场方向变形，并且随着磁场强度的持续增加，液滴的变形将会越来越剧烈。

参考文献

Zhang S T, Niu X D, Li Q P, et al. A numerical investigation on the deformation of ferrofluid droplets[J]. Physics of Fluids, 2023, 35(1): 012102. [百度学术]

Khan A, Niu X D, Li Y, et al. Motion, deformation, and coalescence of ferrofluid droplets subjected to a uniform magnetic field[J]. International Journal for Numerical Methods in Fluids, 2020, 92(11): 1584-1603. [百度学术]

Wu Z G, Troll J, Jeong H H, et al. A swarm of slippery micropropellers penetrates the vitreous body of the eye[J]. Science Advances, 2018, 4(11): eaat4388. [百度学术]

Fan X J, Sun M M, Sun L N, et al. Ferrofluid droplets as liquid microrobots with multiple deformabilities[J]. Advanced Functional Materials, 2020, 30(24): 2000138. [百度学术]

Medina-Sánchez M, Magdanz V, Guix M, et al. Swimming microrobots: soft, reconfigurable, and smart[J]. Advanced Functional Materials, 2018, 28(25): 1707228. [百度学术]

Sun M M, Hao B, Yang S H, et al. Exploiting ferrofluidic wetting for miniature soft machines[J]. Nature Communications, 2022, 13: 7919. [百度学术]

Gao D H, Morley N B, Dhir V. Understanding magnetic field gradient effect from a liquid metal droplet movement[J]. Journal of Fluids Engineering, 2004, 126(1): 120-124. [百度学术]

Korlie M S, Mukherjee A, Nita B G, et al. Modeling bubbles and droplets in magnetic fluids[J]. Journal of Physics Condensed Matter: an Institute of Physics Journal, 2008, 20(20): 204143. [百度学术]

Cunha Lucas H P, Siqueira Ivan R, Oliveira Taygoara F, et al. Field-induced control of ferrofluid emulsion rheology and droplet break-up in shear flows[J]. Physics of Fluids, 2018, 30(12): 122110. [百度学术]

Nochetto R H, Salgado A J, Tomas I. A diffuse interface model for two-phase ferrofluid flows[J]. Computer Methods in Applied Mechanics and Engineering, 2016, 309: 497-531. [百度学术]

Hu Y, Li D C, Niu X D. Phase-field-based lattice Boltzmann model for multiphase ferrofluid flows[J]. Physical Review E, 2018, 98(3): 033301. [百度学术]

Khan A, Niu X D, Li Q Z, et al. Dynamic study of ferrodroplet and bubbles merging in ferrofluid by a simplified multiphase lattice Boltzmann method[J]. Journal of Magnetism and Magnetic Materials, 2020, 495: 165869. [百度学术]

Tang T, Qiao Z H. Efficient numerical methods for phase-field equations[J]. Scientia Sinica Mathematica, 2020, 50(6): 775. [百度学术]

Yue P T, Feng J J, Liu C, et al. A diffuse-interface method for simulating two-phase flows of complex fluids[J]. Journal of Fluid Mechanics, 2004, 515: 293-317. [百度学术]

Dong S, Shen J. A time-stepping scheme involving constant coefficient matrices for phase-field simulations of two-phase incompressible flows with large density ratios[J]. Journal of Computational Physics, 2012, 231(17): 5788-5804. [百度学术]

Jiang G S, Shu C W. Efficient implementation of weighted ENO schemes[J]. Journal of Computational Physics, 1996, 126(1): 202-228. [百度学术]

Huang Z Y, Lin G, Ardekani A M. Consistent, essentially conservative and balanced-force phase-field method to model incompressible two-phase flows[J]. Journal of Computational Physics, 2020, 406: 109192. [百度学术]

Guermond J L, Minev P, Shen J. An overview of projection methods for incompressible flows[J]. Computer Methods in Applied Mechanics and Engineering, 2006, 195(44/45/46/47): 6011-6045. [百度学术]

Rhie C M, Chow W L. Numerical study of the turbulent flow past an airfoil with trailing edge separation[J]. AIAA Journal, 1983, 21(11): 1525-1532. [百度学术]

Huang Z Y, Lin G, Ardekani A M. A mixed upwind/central WENO scheme for incompressible two-phase flows[J]. Journal of Computational Physics, 2019, 387(C): 455-480. [百度学术]

Kim J. A generalized continuous surface tension force formulation for phase-field models for multi-component immiscible fluid flows[J]. Computer Methods in Applied Mechanics and Engineering, 2009, 198(37/38/39/40): 3105-3112. [百度学术]

Liu H H, Valocchi A J, Zhang Y H, et al. Lattice Boltzmann phase-field modeling of thermocapillary flows in a confined microchannel[J]. Journal of Computational Physics, 2014, 256(C): 334-356. [百度学术]

Wang H L, Chai Z H, Shi B C, et al. Comparative study of the lattice Boltzmann models for Allen-Cahn and Cahn-Hilliard equations[J]. Physical Review E, 2016, 94(3): 033304. [百度学术]

Liang H, Shi B C, Guo Z L, et al. Phase-field-based multiple-relaxation-time lattice Boltzmann model for incompressible multiphase flows[J]. Physical Review E, Statistical, Nonlinear, and Soft Matter Physics, 2014, 89(5): 053320. [百度学术]

Xiao Y, Zeng Z, Zhang L Q, et al. A highly accurate bound-preserving phase field method for incompressible two-phase flows[J]. Physics of Fluids, 2022, 34(9): 092103. [百度学术]

Taylor G I. The viscosity of a fluid containing small drops of another fluid[J]. Proceedings of the Royal Society of London Series A, Containing Papers of a Mathematical and Physical Character, 1932, 138(834): 41-48. [百度学术]

Guido S, Villone M. Three-dimensional shape of a drop under simple shear flow[J]. Journal of Rheology, 1998, 42(2): 395-415. [百度学术]

Kennedy M R, Pozrikidis C, Skalak R. Motion and deformation of liquid drops, and the rheology of dilute emulsions in simple shear flow[J]. Computers & Fluids, 1994, 23(2): 251-278. [百度学术]

Cox R G. The deformation of a drop in a general time-dependent fluid flow[J]. Journal of Fluid Mechanics, 1969, 37(3): 601-623. [百度学术]

Hassan M R, Wang C. Magnetic field induced ferrofluid droplet breakup in a simple shear flow at a low Reynolds number[J]. Physics of Fluids, 2019, 31(12): 127104. [百度学术]

Hu Y, Li D C, Jin L C, et al. Hybrid Allen-Cahn-based lattice Boltzmann model for incompressible two-phase flows: the reduction of numerical dispersion[J]. Physical Review E, 2019, 99(2): 023302. [百度学术]

Li Q Z, Lu Z L, Chen Z, et al. An efficient simplified phase-field lattice Boltzmann method for super-large-density-ratio multiphase flow[J]. International Journal of Multiphase Flow, 2023, 160: 104368. [百度学术]

Flament C, Lacis S, Bacri J C, et al. Measurements of ferrofluid surface tension in confined geometry[J]. Physical Review E, 1996, 53(5): 4801-4806. [百度学术]

Home