基于等温伪势格子Boltzmann方法模拟精度提高的研究

Investigation on Improvement of Simulation Accuracy Based on Isothermal Pseudo-Potential Lattice Boltzmann Method

摘要: 采用单组分多相格子Boltzmann方法研究单液滴在等温基板上蒸发时，在壁面边界附近施加虚拟流体层来提高模拟精度是十分必要的。本文采用基于伪势法的流固作用模型，结合构造虚拟流体层的不同方法，研究了不同接触角方案对模拟过程中产生的最大虚假速度和近壁面流体密度偏差的影响，结果表明：基于伪势法构造流固作用力，当液滴接触角θ < 120˚时，采用复制流体密度法构建虚拟流体层会产生最小的虚假速度和密度偏差，且存在密度偏差为0的最佳接触角，是较为理想的构建虚拟流体层，提高模拟精度的接触角方案。

Abstract: It is necessary to apply a ghost fluid layer near the wall boundary to improve the simulation accu-racy of the single-component multiphase lattice Boltzmann method for the evaporation of a single droplet on an isothermal substrate. In this paper, the fluid-solid interaction model based on pseu-do-potential method is used, and the effects of different contact angle schemes on the maximum spurious velocity and the near-solid fluid density deviation are studied in combination with differ-ent methods of constructing the ghost fluid layer. The results show that: The pseudo potential method was used to construct the fluid-solid force, when the droplet contact angle θ < 120˚, the du-plicated fluid density method can produce the smallest spurious velocity and density deviation, and there is an optimal contact angle with zero density deviation. It is an ideal contact angle scheme to construct the ghost fluid layer and improve the simulation accuracy.

1. 引言

润湿现象在自然界和工业中普遍存在，如莲花效应 [1] 、防冻防霜 [2] 和沸腾传热 [3] 等。固体表面的润湿性由接触角来量化，接触角由固体表面的粗糙度和化学均匀性共同决定 [4] 。由于实验技术的局限性，数值模拟对于理解这些复杂过程中润湿性效应的潜在机制起到了非常重要的作用。

格子Boltzmann方法作为一种介观数值方法，具有自动跟踪界面、易于实现复杂几何形状和固体壁面润湿效应的优点。该方法已被广泛应用于研究润湿现象、多相流流动和传热传质问题，如液滴撞击固体表面 [5] 、三相移动接触线 [6] 、多孔介质两相流流动 [7] 和壁面润湿性对沸腾传热的影响 [8] 等。

伪势格子Boltzmann方法(以下简称“LBM”)由于其简单性、通用性和计算效率较高，是应用最广泛的多相LB模型 [9] 。在这种方法中，根据局部流体密度，每个流体节点的“effective mass”模拟了微观流体粒子之间的相互作用。这种方法的优点在于模拟流体具有非理想流体行为，在相变过程中遵循非理想状态方程，并且在相界面处具有表面张力作用，符合真实的物理场。当在伪势多相LBM中的固体壁面边界引入流固相互作用时，可以实现固体壁面的润湿特性 [10] ，固体壁面的接触角由流体粒子间作用力F_cohesive和流体粒子与壁面间的粘附力F_adhesive共同决定，流体粒子间作用力F_cohesive可根据Shan-Chen伪势流体–流体相互作用力计算格式 [9] 或其改进形式 [11] [12] 计算，而流体粒子与壁面间的粘附力F_adhesive大小取决于流体–固体相互作用强度，可根据不同的流固相互作用格式计算 [10] [13] [14] [15] 。

然而在伪势多相LBM中，当对固体壁面施加粘附力F_adhesive时，发现近壁面的流体密度偏离了体系流体密度 [14] [16] [17] [18] 。这种近壁面的流体密度偏差(near-solid fluid density deviation，以下简称“NFDD”)被称为“artificial density distributions” [16] ，或被称为“unphysical mass transfer layer near solid” [18] 。过大的NFDD会对最终模拟结果产生显著的非物理影响。伪势LBM的另一个主要问题是在模拟过程中相界面附近产生非常大的虚假速度(maximum spurious velocity，以下简称“MSV”) [19] ，其会掩盖真实的流体速度，导致算法不稳定。因此，在使用伪势LBM研究等温壁面润湿特性问题时，探究NFDD和MSV的生成影响对于找出最优接触角方案十分重要。

一些研究者 [13] [20] [21] 采用伪势多相LBM实现固体壁面接触角时，建议在壁面附近构造一个附加的虚拟流体层来计算粘附力F_adhesive，但是大多数涉及伪势多相LBM壁面润湿特性的文献都没有明确指出是否在壁面附近添加了的虚拟流体层以及具体的构造方式。此外，在壁面附近添加虚拟流体层对产生MSV和NFDD的影响尚不明确，并且在各个流固相互作用方案中都没有明确提出构建虚拟流体层的最优方法，此方面的研究在文献中也没有得到足够的重视。

因此，本文采用单组分多相伪势LBM，结合三种构造虚拟流体层的方法 [22] ，并与基于伪势的流固作用方案 [15] 结合，以等温壁面液滴蒸发作为媒介，研究了基于伪势流固作用方案中构建虚拟流体层的最优方法，最终确定了产生最小NFDD和MSV的理想接触角方案。

2. 数值方法与问题描述

2.1. 数值方法

在伪势多相格子Boltzmann模型中，每个节点处的流体粒子都具有“effective mass”，x处流体颗粒的“effective mass”可表示为 [9] ：

$\psi \left(\rho \left(x\right)\right)=\sqrt{\frac{2\left(p-\rho \left(x\right)c_s^2\right)}{c_{0}g}}$ (1)

其中 $c_0=6$ ， $g=-1$ ， $\rho \left(x\right)$ 为流体密度， $c_s=1/\sqrt{3}$ ，压力P可由下式P-R状态方程计算得出 [11] ：

$p=\frac{\rho RT}{1-b\rho }-\frac{a{\rho }^2\epsilon \left(T\right)}{1+2b\rho -b^2{\rho }^2}$ (2a)

其中 $a=0.4572R^2{T}_c^2/p_c$ $b=0.0778R{T}_c/p_c$ ，且 $\epsilon \left(T\right)$ 可通过下式计算：

$\epsilon \left(T\right)={\left[1+\left(0.3746+1.5423\omega -0.2699{\omega }^2\right)\left(1-\sqrt{T/T_c}\right)\right]}^2$ (2b)

其中Tc和Pc表示临界温度和临界压力，本文中水的ω为0.344，F_cohesive可由下式计算得到：

$F_{\text{cohesive}}\left(x\right)=-\beta \psi \left(x\right){\displaystyle \underset{x^{\prime }}{\sum }G\left(x,x^{\prime }\right)\psi \left(x^{\prime }\right)}\left(x^{\prime }-x\right)-\frac{1-\beta }{2}{\displaystyle \underset{x^{\prime }}{\sum }G\left(x,x^{\prime }\right){\psi }^2\left(x^{\prime }\right)\left(x^{\prime }-x\right)}$ (3)

方程(3)中的 $G\left(x,x^{\prime }\right)$ 为格林函数，满足 $G\left(x,x^{\prime }\right)=G\left(x^{\prime },x\right)$ ，且 $\beta =1.16$ 。

格子Boltzmann方法利用流体粒子微团的分布函数 $f_i\left(x,t\right)$ 来表征计算区域内流体的运动，具有LBGK碰撞算子的 $f_i\left(x,t\right)$ 的单弛豫控制方程表示为：

$f_i\left(x+e_i\Delta t,t+\Delta t\right)-f_i\left(x,t\right)=-\frac{1}{\tau }\left[f_i\left(x,t\right)-f_i^{eq}\left(x,t\right)\right]+\Delta f_i\left(x,t\right)$ (4)

其中 $f_i\left(x,t\right)$ 为𝑡时刻位于𝑥位置的速度为 $e_i$ 的流体粒子微团的分布函数， $\Delta t$ 是时间步长，一般取 $\Delta t=1$ ，t为与流体运动黏度有关的无量纲松弛时间， $f_i^{eq}\left(x,t\right)$ 是流体粒子微团平衡态分布函数， $\Delta f_i\left(x,t\right)$ 为体积力项，由下式计算得到：

$\Delta f_i\left(x,t\right)=f_i^{eq}\left[\rho \left(x,t\right),u\left(x,t\right)+\Delta u\right]-f_i^{eq}\left[\rho \left(x,t\right),u\left(x,t\right)\right]$ (5)

其中 $\Delta u=F\left(x\right)\cdot \Delta t/\rho$ ， $F\left(x\right)$ 为合力，可表示为：

$F\left(x\right)=F_{\text{cohesive}}\left(x\right)+F_{\text{adhesive}}\left(x\right)$ (6)

其中F_adhesive表示的是流体与壁面之间的粘附力，其在数值上等于流固相互作用力F_s与流体–虚拟流体相互作用力F_ghost之和，即：F_adhesive = F_s + F_ghost。

流固相互作用力F_s可根据以下流固相互作用模型计算 [15] ：

$F_{\text{s}}\left(x\right)=-G_s\psi \left(x\right){\displaystyle \underset{i}{\sum }{\omega }_{i}s\left(x+e_i\right)e_i}$ (7)

流体–虚拟流体相互作用力F_ghost可通过以下三种方式构建 [22] ：

A. 无虚拟流体层：

$F_{\text{ghost}}\left(x\right)=0$ (8)

B. 复制法构建虚拟流体层：

$F_{\text{ghost}}\left(x\right)=-\beta \psi \left(x\right){\displaystyle \underset{i}{\sum }G\left(x,x+e_i\right)\psi \left(x_{\text{neighbor}}\right)e_i}-\left[\left(1-\beta \right)/\text{2}\right]{\displaystyle \underset{i}{\sum }G\left(x,x+e_i\right){\psi }^2\left(x_{\text{neighbor}}\right)e_i}$ (9)

其中对于壁面， $\left(x_{\text{neighbor}},y_{\text{neighbor}}\right)=\left(x_{\text{solid}},y_{\text{solid}}+1\right)$ 。

C. 常数法构建虚拟流体层：

$F_{\text{ghost}}\left(x\right)=-\beta \psi \left(x\right){\displaystyle \underset{i}{\sum }G\left(x,x+e_i\right)\psi \left({\rho }_{\mathrm{int}}\right)e_i}-\left[\left(1-\beta \right)/\text{2}\right]{\displaystyle \underset{i}{\sum }G\left(x,x+e_i\right){\psi }^2\left({\rho }_{\mathrm{int}}\right)e_i}$ (10)

其中 ${\rho }_{\mathrm{int}}=\left({\rho }_{\text{liquid,sat}}+{\rho }_{\text{vapor,sat}}\right)/\text{2}$ 。

D. 当地流体密度法构建虚拟流体层：

$F_{\text{ghost}}\left(x\right)=-\beta \psi \left(x\right){\displaystyle \underset{i}{\sum }G\left(x,x+e_i\right)\psi \left(x\right)e_i}-\left[\left(1-\beta \right)/\text{2}\right]{\displaystyle \underset{i}{\sum }G\left(x,x+e_i\right){\psi }^2\left(x\right)e_i}$ (11)

Table 1. Contact angle schemes in this paper

表1. 本文采用的接触角方案

2.2. 问题描述

图1给出了采用表1总结的接触角方案，模拟了等温基板的液滴蒸发过程，具体参数设置如图所示。根据Maxwell构造的汽液共存曲线 [15] ，可以得到饱和液相密度ρ_l = 5.91，饱和汽相密度ρ_v= 0.58。饱和液相与汽相的运动粘度分别为υ_l = 0.17，υ_v = 1，液滴的初始半径R = 50，等温基板厚度H = 20∆y，计算区域为 $NX\times NY$ ，经网格独立性验证后，计算区域选择为400 × 150，计算域上下边界为无滑移边界条件，左右边界为周期性 $NX\times NY$ 边界条件，固体壁面设置为反弹边界条件。

Figure 1. Diagram of droplet evaporation on isothermal substrate

图1. 等温基板液滴蒸发示意图

3. 模型验证

为了验证所使用的LBM模型的正确性，采用拉普拉斯定律对模型进行了可靠性检验。本文将单个液滴置于计算域的中心，边界条件设置为周期性边界条件，模拟所选取的液滴半径为20、25、30、35、40、45、50格子单位，结果如图2所示。从图中可以看出，当表面张力一定时，液滴内外压差与1/R呈线性关系，符合Laplace定律。

本文通过对液滴蒸发的D²定律 [23] 的分析，进一步验证了模型的准确性。在模拟过程中将液滴放在计算区域中心，采用对流与周期性边界条件，液滴温度设置为饱和温度，其中D₀为液滴初始直径，D为液滴蒸发后的直径，图3可以看出 ${\left(D/D_0\right)}^2$ 与蒸发时间呈线性关系，符合已知的D²定律，且与Nishiwaki [24] 的实验结果相吻合。

Figure 2. Verification of Laplace law

图2. Laplace定律验证图像

Figure 3. Verification of D² law

图3. D²定律验证图像

4. 结果与讨论

采用基于伪势法的流固作用格式结合四种虚拟流体层的构建方法(表1中的方法1A，1B，1C和1D)，对等温基板固着液滴蒸发过程中产生的MSV和NFDD等问题进行探究，在每一个接触角模型中，流固作用力F_s由方程(7)给出，F_ghost有方程(8~11)给出，最终相对理想的接触角方案将在方案1A，1B，1C和1D中确定。

4.1. 最大虚假速度(MSV)的研究

图4表明基于伪势法–不施加虚拟流体层(1A)和基于伪势法–常流体密度法(1C)，在壁面不同接触角

下两种方案产生的最大虚假速度(MSV， ${\left|u_{\text{sp}}\right|}_{\max }$ )十分相近，且相较于其他虚拟流体层的构建方案，保持

在最小数值；以复制流体密度法(1B)构建虚拟流体层时，在接触角θ < 120˚时MSV维持在较小数值，但是当接触角θ > 120˚时，MSV会急剧增大，可能最终会导致算法的不稳定；以当地流体密度法(1D)构建虚拟流体层时，此方案的MSV相较于其他方案存在着显著提高，并且MSV会随着接触角的增加而变大。综上所述，从图4可以得出，在基于伪势的流固作用方案中，无虚拟流体层(A)、复制法(B)和常数法(C)构造虚拟流体层具有最小的MSV，因此在接触角θ < 120˚，在等温固体壁面上使用以上三种方法构造虚拟流体层是数值最稳定的模拟方法。

Figure 4. Comparison of MSV in different ghost fluid layer construction methods in the pseudo-potential-based interaction scheme at different contact angles θ (1A: no-ghost fluid layer method; 1B: duplicated method; 1C: constant value method; 1D: local fluid density method)

图4. 基于伪势流固作用方案在不同接触角θ下以不同虚流体层构造方法的MSV比较(1A：无虚流体层；1B：复制法；1C：常数法；1D：当地流体密度法)

4.2. 液相/汽相近壁面密度分布的研究

图5(a)，图5(b)，图5(c)，图5(d)分别给出了采用不同方案构造虚拟流体层(1A，1B，1C和1D)，在不同接触角θ下，从体系流体到近壁面(y = 21~30)处的液相流体密度(黑色实线)和汽相流体密度(黑色虚线)的分布情况，其中液相或汽相的体系流体密度在图中用红色水平点划线表示。

当等温壁面无虚拟流体层时(图5(a)中的方案1A)，固体壁面上的近壁面液相密度相较于体系液相密度偏离的流体层厚度为6∆y (y = 21~26)，近壁面汽相密度相较于体系汽相密度偏离的流体层厚度为8∆y (y = 21~28)，并且随着接触角的增大，近壁面液相流体密度偏离值会逐渐变大，汽相流体密度偏离值会逐渐缩小；同时，当等温壁面采用常流体密度法构建虚拟流体层时(图5(c)中的方案1C)，其近壁面产生的流体密度相较于体系流体的偏离程度的变化规律与方案1A相似；当等温壁面采用复制流体密度法构建虚拟流体层时(图5(b)中的方案1B)，此时近壁面液相密度相较于体系液相密度偏离的流体层厚度为4∆y (y = 21~24)，汽相密度相较于体系汽相密度偏离的流体层厚度为7∆y (y = 21~27)，可以得出结论，存在着密度偏离度为0的最佳接触角θ = 90˚，当接触角偏离θ = 90˚时，近壁面液相/汽相流体的密度偏离程度会突然增大；同时，当等温壁面采用当地流体密度法构建虚拟流体层时(图5(d)中的方案1D)，近壁面近壁面液相/汽相流体的密度偏离程度变化有着与方案1B相似的规律，且同样存在着密度偏离度为0的最佳接触角θ = 65.1˚。

值得注意的是，当采用方案1B和方案1C构造虚拟流体层时，密度偏离度为0对所应的接触角大小虽然不同，但这两种方案所对应的调节接触角的流固作用系数相同：G_s = 0，可以得出结论，采用基于伪势法的流固作用格式，以复制流体密度法(1B)或当地流体密度法(1D)构造虚拟流体层，当流固作用系数G_s = 0时，可以将近壁面的密度偏离度完全消除。

(a) 接触角方案1A条件下近壁面液相(汽相)密度分布

(b) 接触角方案1B条件下近壁面液相(汽相)密度分布

(c) 接触角方案1C条件下近壁面液相(汽相)密度分布(d) 接触角方案1D条件下近壁面液相(汽相)密度分布

Figure 5. Effects of different ghost fluid layer construction methods on near-wall fluid density distribution using the pseudo-potential-based interaction scheme: (a) no-ghost fluid layer method (1A), (b) duplicated method (1B), (c) onstant value method (1C), (d) local fluid density method (1D)

图5. 基于伪势流固作用格式的不同虚流体层构造方法对近壁流体密度分布的影响：(a) 无虚流体层(1A)，(b) 复制法(1B)，(c) 常数法(1C)，(d) 当地流体密度法(1D)

4.3. 近壁面密度偏差(NFDD)的研究

当采用不同的接触角方案(1A，1B，1C和1D)实现等温壁面的润湿特性时，如图5所示，近壁面的液相/汽相的流体层密度会相较于体系流体密度发生非物理性偏离，为了定量评估这种非物理性流体密度偏差，将近壁面液相(蒸相)的密度偏差因子定义为：

${\epsilon }_{\text{liquid}}=\left|\rho \left(NX/\text{2},H+1\right)-\rho \left(NX/\text{2},H+10\right)\right|/\rho \left(NX/\text{2},H+10\right)$ (12)

${\epsilon }_{\text{vapor}}=\left|\rho \left(0.84\times NX,H+1\right)-\rho \left(0.84\times NX,H+10\right)\right|/\rho \left(0.84\times NX,H+10\right)$ (13)

其中 $\rho \left(NX/\text{2},H+1\right)$ 表示的是x = NX/2处平壁附近流体层(y = H + 1)的液相密度， $\rho \left(NX/\text{2},H+10\right)$ 表示的是x = NX/2处体系液相流体密度，此时体系液相流体密度选取在y = H + 10，这是因为从图5可以看出，流体的密度偏差只发生在近壁面附近的流体层内(流体层厚度 < 10∆y)，选择近壁面汽相密度值 $\rho \left(0.84\times NX,H+1\right)$ 和体系汽相流体密度值 $\rho \left(0.84\times NX,H+10\right)$ 的原因也与液相流体密度值的选择原因类似。

图6表明采用方案1A和1C，近壁面液相密度偏差 ${\epsilon }_{\text{liquid}}$ 随着接触角的增加而明显增大，汽相密度偏差 ${\epsilon }_{\text{vapor}}$ 随接触角的增大而减小，所以在任意接触角下，采用基于伪势的流固作用格式时，不推荐使用无虚拟流体层(A)和常流体密度法(C)来构建虚拟流体层；与之相反，采用复制流体密度法(B)和常流体密度法(D)构建流体层， ${\epsilon }_{\text{liquid}}$ 和 ${\epsilon }_{\text{vapor}}$ 相较于方案1A和1C要小得多。通过比较图5与图6得出的结果，可以得出结论：采用基于伪势法的流固作用格式，在接触角θ < 120˚时，以复制流体密度法(B)构造虚拟流体层是最合理的方案，且存在最佳接触角θ = 90˚。

Figure 6. Comparison of near-wall liquid density deviation ${\epsilon }_{\text{liquid}}$ and near-solid vapor density deviation ${\epsilon }_{\text{vapor}}$ in different ghost fluid layer construction method by using the pseudo-potential-based interaction scheme at different contact angles θ (1A: no-ghost fluid layer method; 1B: duplicated method; 1C: constant value method; 1D: local fluid density method)

图6. 基于伪势法不同构造虚拟流体层方式的 ${\epsilon }_{\text{liquid}}$ 和 ${\epsilon }_{\text{vapor}}$ 随接触角θ的变化情况(1A：无虚流体层；1B：复制法；1C：常数法；1D：当地流体密度法)

5. 结论

本文采用单组分多相格子Boltzmann方法，主要研究了以不同方式施加虚拟流体层对等温基板液滴蒸发模拟精度的影响，主要结论如下：

1) 基于伪势法构造流固作力，以复制流体密度法(B)施加虚拟流体层，在接触角θ < 120˚时产生的MSV保持在较小数值，但不适合模拟较疏水的润湿表面(θ > 120˚)，会使MSV急剧增大，而以常流体密度法和当地流体密度法施加虚拟流体层，产生的MSV始终维持在较小的数值。

2) 采用基于伪势法的流固作用格式，在接触角θ < 120˚时，以复制流体密度法(B)构造虚拟流体层产生的NFDD最小，且存在壁面密度偏差为0的最佳接触角，是较为理想的构建虚拟流体层，提高模拟精度的接触角方案。

NOTES

^*通讯作者。