变温场诱导相分离的格子Boltzmann数值模拟

Lattice Boltzmann Numerical Simulation of Phase Separation Induced by Variable Temperature Field

摘要: 相分离是环境条件发生变化导致原来混合流体分离出两相或多相的不稳定倾向和过程，是多相流研究的核心内容之一。本文采用自由能格子Boltzmann模型，将流体黏度与温度变化过程耦合，研究二元混合流体在温度变化条件下的相分离行为。

Abstract: Phase separation is an unstable tendency and process of separating two or more phases from the original mixed fluids due to the change in environmental conditions. It is one of the core contents in the study of multiphase fluids. In this paper, the phase separation behavior of binary fluids under the change of temperature is studied by the free-energy lattice Boltzmann model coupling the fluid viscosity and temperature.

1. 引言

变温场相分离是二元流体许多理论和工业生产研究的重要课题 [1] [2] [3] [4] [5]。在很多科学研究和实际工程中，要求以受控方式进行相分离过程，以创建规则结构。例如，在二元合金 [6] 中，慢冷却用于产生不同材料交替带的最佳序列。在聚合物混合物中，分层形态应通过适当的热驱动进行控制 [7] [8]。目前，变温场诱导相分离的研究仍然十分活跃，并且应用于材料、化学以及生物学等领域 [9] [10] [11] [12]。一些学者已经对相分离问题进行了数值模拟研究，包括使用有限差分法、有限体积法等传统数值方法。近年来，随着格子Boltzmann方法(LBM)的发展，该方法也被应用于研究相分离问题。格子Boltzmann方法是一种介观方法，能够刻画更多的物理细节。此外，与传统数值方法相比，LBM具有数值优势，例如算法设计简单，易于编程以及能够处理复杂的边界等。本文建立了一个耦合流体黏度与温度变化过程的格子Boltzmann模型来模拟二元混合流体在冷边界条件下的相分离，并进一步阐明相分离的形成与发展机理。

2. 相离分数学模型

2.1. 耦合温度与粘度的多相模型

对于二元混合流体，以总密度 $\rho$ 、浓度差 $\phi$ 表示的质量和动量守恒方程 [13]，和宏观对流扩散方程，其形式如下

$\frac{\partial \rho }{\partial t}+\nabla \left(\rho \stackrel{\to }{u}\right)=0$ (1)

$\frac{\partial \left(\rho \stackrel{\to }{u}\right)}{\partial t}+\nabla \left(\rho \stackrel{\to }{u}\stackrel{\to }{u}\right)=-\nabla P_{\alpha \beta }+{\eta }_{mix}{\nabla }^2\left(\rho \stackrel{\to }{u}\right)+\nabla \left(\lambda \left(\rho \right)\nabla \left(\rho \stackrel{\to }{u}\right)\right)$ (2)

$\frac{\partial \varphi }{\partial t}+\nabla \left(\phi \stackrel{\to }{u}\right)=\Gamma \theta {\nabla }^2\delta \mu -\theta \nabla \left(\frac{\phi }{\rho }\nabla P_{\alpha \beta }\right)$ (3)

$\frac{\partial T}{\partial \tilde{t}}+\nabla \left(T\stackrel{\to }{u}\right)=\nabla \cdot \left(D_{mix}\nabla T\right)+F$ (4)

式中， $P_{\alpha \beta }$ 、 $\stackrel{\to }{u}$ 、 ${\eta }_{mix}$ 、 $\lambda \left(\rho \right)$ 、 $\Gamma \theta$ 、 $D_{mix}$ 和F分别是二元流体中混合物的压力张量、速度、混合运动粘度、体粘度、宏观迁移率、热扩散率和外力 [14]。结果表明，参数 $\Gamma$ 是与两相流体流动性相关的系数， $\theta$ 是宏观输运系数 [15] [16]。本文将宏观迁移率参数 $\Gamma \theta$ 设为0.2，以获得数值稳定性。混合粘度 ${\eta }_{mix}$ 和热扩散系数 $D_{mix}$ 可以写为

${\eta }_{mix}={\eta }_A\cdot \frac{\rho +\phi }{2\rho }+{\eta }_B\cdot \frac{\rho -\phi }{2\rho }$ (5)

$D_{mix}=D_A\cdot \frac{\rho +\phi }{2\rho }+D_B\cdot \frac{\rho -\phi }{2\rho }$ (6)

其中， $D_A$ 、 $D_B$ 和 ${\eta }_A$ 、 ${\eta }_B$ 分别是组分A和B的热扩散系数和初始粘度。二元体系的温度依赖粘度 $\eta$ [17] 由下式确定

$\eta =\frac{{\eta }_{mix}}{1+\beta \left(T-T_{\infty }\right)}$ (7)

其中 $T_{\infty }$ 是初始温度。温度粘度系数 $\beta$ 是恒定的，取决于流体的热性质，它表示与温度相关的粘度灵敏度。这里值得一提的是，当 $\beta \to 0$， $\eta \to {\eta }_{mix}$ (常数)。还需要注意的是， $\beta$ 对于液体为正，对于气体为负。宏观方程式(2)修改为

$\frac{\partial \left(\rho \stackrel{\to }{u}\right)}{\partial t}+\nabla \left(\rho \stackrel{\to }{u}\stackrel{\to }{u}\right)=-\nabla P_{\alpha \beta }+{\eta }_{mix}{\nabla }^2\left(\rho \stackrel{\to }{u}\right)+\nabla \left(\lambda \left(\rho \right)\nabla \left(\rho \stackrel{\to }{u}\right)\right)-\nabla \left(\frac{{\eta }_{mix}\beta \left(T-T_{\infty }\right)}{1+\beta \left(T-T_{\infty }\right)}\nabla \left(\rho \stackrel{\to }{u}\right)\right)$ (8)

这些方程描述了非均匀温度系统的动力学。

通过压力张量 $P_{\alpha \beta }$ 和化学势差 $\delta \mu$ [14] [15]，混合物的热力学由自由能确定：

$\psi ={\displaystyle \int \left(\varphi \left(T,\rho ,\phi \right)+\frac{\kappa }{2}{\left(\nabla \rho \right)}^2+\frac{\kappa }{2}{\left(\nabla \phi \right)}^2\right)}$ (9)

温度 $T$ 下的体自由能密度为

$\varphi \left(T,\rho ,\phi \right)=\frac{\lambda }{4}\rho \left(1-\frac{{\phi }^2}{{\rho }^2}\right)-T\rho +\frac{T}{2}\left(\rho +\phi \right)\ln \left(\frac{\rho +\phi }{2}\right)+\frac{T}{2}\left(\rho -\phi \right)\ln \left(\frac{\rho -\phi }{2}\right)$ (10)

其中 $\lambda$ 是二元液体的相互作用强度， $\kappa$ 是界面宽度。当温度T降至 $T_c=\frac{\lambda }{2}$ 时，混合流体被分为两个相，密度差为 $\pm \phi$。化学势差和压力张量 [18] 为

$\delta \mu \left(\phi ,\rho ,T\right)=-\frac{\lambda }{2}\frac{\phi }{\rho }+\frac{T}{2}\ln \left(\frac{\rho +\phi }{\rho -\phi }\right)-\kappa {\nabla }^2\phi -\kappa T\left(\nabla \phi \nabla \frac{1}{T}\right)$ (11)

和

$P_{\alpha \beta }=p_0{\delta }_{\alpha \beta }+\kappa \left(\frac{\partial \rho }{\partial x_{\alpha }}\frac{\partial \rho }{\partial x_{\beta }}+\frac{\partial \phi }{\partial x_{\alpha }}\frac{\partial \phi }{\partial x_{\beta }}\right)$ (12)

其中

$p_0=\rho T-\frac{\kappa }{2}\left(\rho {\nabla }^2\rho +\phi {\nabla }^2\phi \right)-\frac{\kappa }{2}\left({\left|\nabla \rho \right|}^2+{\left|\nabla \phi \right|}^2\right)$

在最常见的Raylareigh-Benard对流形式中，粘性流体被限制在两个保持不同温度场的水平刚性边界之间。当流体具有位置热膨胀系数，且重力与温度梯度方向相同时，净浮力与重力方向相反 [19]。当两个边界之间的温差超过某个阈值时，静态导热状态变得不稳定，对流突然发生。众所周知的Boussinesq近似通常用于自然对流的研究。在这种近似下，假设材料特性与温度无关，但体积力项除外，其中流体密

度 $\rho$ 被假设为温度的线性函数，即 $\frac{\rho }{{\rho }_{\infty }}=1+D\left(T-T_{\infty }\right)$。这里， ${\rho }_{\infty }$ 和 $T_{\infty }$ 分别为参考点处的密度和温度，

D为恒定热膨胀系数。因此，重力为 $F={\rho }_{\infty }g+{\rho }_{\infty }gD\left(T-T_{\infty }\right)$。在将第一项吸收到压力中后，有效体积力与温度变化成线性比例。在我们的模拟中，当压缩可以忽略不计时(如在小马赫数限制下)，速度场近似为扩散力，温度场满足被动标量方程(4)。

2.2. 二维LBGK模型

我们使用二维九速度(D2Q9)模型将Navier-Stokes方程(NS)与二维(2D) Boltzmann格子Bhatnagar- Gross-Krook模型(LBGK)对应 [20]。网格上的粒子速度 ${\stackrel{\to }{e}}_i$ 如图1所示，可以写为：

${\stackrel{\to }{e}}_i=\{\begin{array}{l}\left(0,0\right)i=0\\ c\left(\cos \left(i-1\right)\frac{\pi }{2},\sin \left(i-1\right)\frac{\pi }{2}\right)i=1,\cdots ,4\\ \sqrt{2}c\left(\cos \left(2i-1\right)\frac{\pi }{4},\sin \left(2i-1\right)\frac{\pi }{4}\right)i=5,\cdots ,8\end{array}$ (13)

其中 $c=\frac{\delta x}{\delta t}$，这里 $\delta x$ 为空间步长(格子边长)， $\delta t$ 为时间步长。

Figure 1. Schematic diagram of D2Q9 velocity model

图1. D2Q9速度模型示意图

通过LBM，我们使用双分布函数模型来描述二元流体中的相分离和温度场。分别通过两组分布函数 $f_i$ 和 $g_i$，在空间和时间上离散二元流体混合系统的密度场和相分离场的演化方程，每组份流体与速度矢量 ${\stackrel{\to }{e}}_i$ 相关：

$f_i\left(\stackrel{\to }{x}+{\stackrel{\to }{e}}_i\delta t,t+\delta t\right)-f_i\left(\stackrel{\to }{x},t\right)=-\frac{1}{{\tau }_f}\left(f_i\left(\stackrel{\to }{x},t\right)-f_i^{eq}\left(\stackrel{\to }{x},t\right)\right)$ (14)

$g_i\left(\stackrel{\to }{x}+{\stackrel{\to }{e}}_i\delta t,t+\delta t\right)-g_i\left(\stackrel{\to }{x},t\right)=-\frac{1}{{\tau }_g}\left(g_i\left(\stackrel{\to }{x},t\right)-g_i^{eq}\left(\stackrel{\to }{x},t\right)\right)$ (15)

其中 $\stackrel{\to }{x}$ 为网格节点位置， ${\tau }_f$ 和 ${\tau }_g$ 是松弛时间， $f_i^{eq}$ 和 $g_i^{eq}$ 是平衡态分布函数。

通过如下的粒子分布函数，现有的热LBM可以模拟温度场的演变，

$h_i\left(\stackrel{\to }{x}+{\tilde{e}}_i\delta \tilde{t},t+\delta \tilde{t}\right)-h_i\left(\stackrel{\to }{x},t\right)=-\frac{1}{{\tau }_h}\left(h_i\left(\stackrel{\to }{x},t\right)-h_i^{eq}\left(\stackrel{\to }{x},t\right)\right)$ (16)

其中 ${\tau }_h$ 是与热扩散系数 $D_A=D_B=\left({\tau }_h-0.5\right){\tilde{c}}^2\delta \tilde{t}/3$ 相关的无量纲松弛时间。 $\delta \tilde{t}$ 是分布函数 $h_i$ 的时间步长，且 $\tilde{c}=\frac{\delta x}{\delta \tilde{t}}$，

${\tilde{e}}_i=\{\begin{array}{l}\left(0,0\right)i=0\\ \tilde{c}\left(\cos \left(i-1\right)\frac{\pi }{2},\sin \left(i-1\right)\frac{\pi }{2}\right)i=1,\cdots ,4\\ \sqrt{2}\tilde{c}\left(\cos \left(2i-1\right)\frac{\pi }{4},\sin \left(2i-1\right)\frac{\pi }{4}\right)i=5,\cdots ,8\end{array}$ (17)

宏观流速 $\stackrel{\to }{u}$，密度 $\rho$，密度差 $\phi$ 和温度T关于平衡态分布函数 $f_i^{eq}$， $g_i^{eq}$ 和 $h_i^{eq}$ 的表达式如下：

$\rho ={\displaystyle \underset{i}{\sum }{f}_i^{eq}}$ (18)

$\rho u_{\alpha }={\displaystyle \underset{i}{\sum }{f}_i^{eq}}{e}_{i\alpha }$ (19)

$\phi ={\displaystyle \underset{i}{\sum }{g}_i^{eq}}$ (20)

$T={\displaystyle \underset{i}{\sum }{h}_i}$ (21)

平衡态密度分布函数的高阶矩定义如下：

${\displaystyle \underset{i}{\sum }{f}_i^{eq}}{e}_{i\alpha }{e}_{i\beta }=P_{\alpha \beta }+\rho u_{\alpha }{u}_{\beta }$ (22)

${\displaystyle \underset{i}{\sum }{g}_i^{eq}}{e}_{i\alpha }=\phi u_{\alpha }$ (23)

${\displaystyle \underset{i}{\sum }{g}_i^{eq}}{e}_{i\alpha }{e}_{i\beta }=\Gamma \delta \mu {\delta }_{\alpha \beta }+\phi u_{\alpha }{u}_{\beta }$ (24)

${\displaystyle \underset{i}{\sum }{h}_i^{eq}}{e}_{i\alpha }=T{u}_{\alpha }$ (25)

这些矩必须满足宏观方程(1)~(3)，以描述二元液体混合物的动力学。这里，自由能格子Boltzmann格式也适用于方程组(22)~(24)。通过Chapman Enskog展开，可以通过演化方程(14)~(15)恢复连续方程、NS方程和相场方程(1)、(3)和(8)，从而得到

$\theta =c^2\delta t\left({\tau }_f-\frac{1}{2}\right)$

${\eta }_{mix}=\frac{\left(2{\tau }_f-1\right)\left(1+\beta \left(T-T_{\infty }\right)\right)}{6}{c}^2\delta t$

$\lambda \left(\rho \right)=\left({\tau }_f-\frac{1}{2}\right)\left(\frac{2c^2}{3}-\frac{\text{d}{p}_0}{\text{d}\rho }\right)\delta t$ (26)

来自等式(14)和(15)的这一推导，受参考文献 [16] 的启发。有必要根据计算的温度依赖粘度 $\frac{{\eta }_{mix}}{1+\beta \left(T-T_{\infty }\right)}$ 在每个网格点更新 ${\tau }_f$。但需要注意，对于 $\beta$ 在0~9.09范围内都保证LBGK的稳定性。

2.3. 边界处理

在所有模拟中，在上边界和下边界施加了非滑移边界条件，正如针对二元剪切流体系统提出的那样。该算法可以使研究人员严格保持边界墙上的质量和动量。可用于模拟温度场的热边界条件 [21] 表示为局部已知值和边界处的修正值的组合：

$h_i\left(\stackrel{\to }{x},t\right)=h_i^{\ast }\left(\stackrel{\to }{x},t\right)+{\omega }_i{G}_c$ (27)

其中 $G_c$ 是内能影响的修正值。Dirichlet热边界条件适用于温度T恒定的边界节点。对于迁移后的上边界，在上排的每个位置都知道函数 $h_0$ 、 $h_1$ 、 $h_3$ 、 $h_2$ 、 $h_5$ 和 $h_6$。然后使用方程(27)确定 $h_4$ 、 $h_7$ 和 $h_8$，要求上边界温度为 $T_c$，这也适用于下边界。修正后上边界分布函数为

$G_c=\frac{T_c-\left(h_0+h_1+h_3+h_2+h_5+h_6\right)}{{\omega }_4+{\omega }_7+{\omega }_8}$ (28)

$h_4=h_4^{\ast }+{\omega }_4{G}_c$

$h_7=h_7^{\ast }+{\omega }_7{G}_c$

$h_8=h_8^{\ast }+{\omega }_8{G}_c$ (29)

其中分布函数 $h_4^{\ast }$， $h_7^{\ast }$ 和 $h_8^{\ast }$ 由上一步的计算值决定，即 $h_i^{\ast }\left(\stackrel{\to }{x},t\right)=h_i\left(\stackrel{\to }{x},t-\delta \tilde{t}\right)$。

3. 数值模拟结果与分析

我们演示耦合温度与粘度的自由能格子Boltzmann模型在 $\rho =1.0$， $\phi =0$ 和 $\stackrel{\to }{u}=0$ 情况下的简单应用。图2给出了混合流体初始状态示意图。在模拟中， $\rho$ 的值在系统的整个时间演化过程中保持不变。 $\phi$ 的初始值对应于两种流体完全混合的对称系统。为便于描述热扩散和对流之间的竞争，引入Prandtl数(Pr)，

Figure 2. Initial state diagram of mixed fluid

图2. 混合流体初始状态示意图

其定义为 $\mathrm{Pr}={\eta }_{\infty }/D$。其它计算参数为 $k=6.67\times {10}^{-3}$， ${\tau }_g=\left(1+1/\sqrt{3}\right)/2$， $\delta x=1.0$， $\Gamma \theta =0.2$。格子数为 $256\times 256$，初始温度为 $T_{\infty }=0.6$ 和边界温度 $T_C=0.49$ (对应于最大密度差 ${\phi }_C\approx \pm 0.55$ )，参数 $\beta$ 与式(7)

相关，其取值范围为 $\left[0,\frac{1}{T_{\infty }-T_C}\right]$，这里 $\beta ={10}^{-2}$，热扩散系数 $D_A=D_B=0.1$，粘度 ${\eta }_A={\eta }_B=0.1$ ( $\mathrm{Pr}=1$ )。

图3给出了在温度变化条件下的相分离过程。从图3(a)可以看到，在上下低温边界附近首先出现相分离，这是因为在靠近低温边界地方通过传热使得温度首先降低，从而导致原来混合的流体发生分离。随着时间演化，内部温度高的地方不断向低温边界方向传热，相分离区域不断向内部推进(图3(b)~(e))。在这个过程中，由于温差引起的对流也对传热和相分离产生影响，加快传热和相分离的进程。此外，经过一段时间后，表面张力引起的速度场可以加速同相聚集的过程。经过较长的一段时间后，相分离渐渐

(a) t = 5000 s (b) t = 15,000 s (c) t = 30,000 s (d) t = 45,000 s (e) t = 60,000 s (f) t = 75,000 s (g) t = 100,000 s (h) t = 150,000 s

Figure 3. Evolution of phase separation

图3. 相分离的演化过程

趋于稳态，聚集和破碎两种机制将达到平衡(图3(f)~(h))。图3相分离的数值模拟结果与文献 [22] 的结果定性一致。在这个相分离过程中，经历两个阶段，一个阶段是亚稳相分解阶段，另一个阶段是相畴增长阶段。亚稳相分解阶段是指发生相分离的前阶段，在这个阶段，相分离已经发生，但还没有形成明显的两相聚集。相畴增长阶段是指发生相分离的后阶段，在这个阶段，相分离产生了明显的两相聚集，在温度变化和表面张力的共同作用下，两相聚集的液滴不断长大。

4. 总结

本工作采用耦合温度和黏度的自由能格子Boltzmann模型研究了由温度变化引起的二元混合流体的相分离过程。在低温边界附近先发生相分离，随着传热的进行，相分离渐渐向内部高温区域发展，最后相分离趋于稳定。相分离经历亚稳相分解和相畴增长两个阶段。计算结果还表明，耦合温度和黏度的自由能格子Boltzmann模型能够有效模拟变温场诱导相分离过程。

基金项目

国家自然科学基金(11804355, 31800083)。

参考文献