固壁上液滴动力学的LBM建模与计算

LBM Modeling and Calculation for Droplet Dynamics on Solid Wall

摘要: 本文采用Shan-Chen多相格子Boltzmann方法建立液滴在固壁上的动力学模型。并根据该模型模拟了液滴在水平固壁上的铺展行为以及在垂直固壁上的滑移行为。

Abstract: In this paper, the dynamic model for the droplet on the solid wall is studied with Shan-Chen multiphase lattice Boltzmann method. Based on the model, the spreading behavior of droplet on the horizontal solid wall and the slip behavior on the vertical solid wall are simulated.

1. 引言

固液相互作用流在许多领域都有应用，如地质封存二氧化碳、地下污染物修复、石油开采和质子交换膜燃料电池(PMEFC)等。

当面对这些问题时，数值技术可以提供很大帮助。典型的方法是基于离散化Navier-Stokes方程，但不容易模拟细尺度问题，且非线性对流项的处理存在困难。

格子Boltzmann方法(LBM)是一种介观尺度方法，它能够正确处理微尺度的复杂性，同时再现流动的宏观特征。Shan-Chen模型是一种多相LBM [1] ，已成功应用于许多研究 [2] 。近年来，流固相互作用的研究引起广泛关注 [3] [4] [5] [6] [7] 。

本文根据LBM的介观特性，利用Shan-Chen多相模型建立液滴在固壁上的动力学模型。本文的主要内容如下：首先对LBM和SC模型进行了描述，接着建立固壁与液滴相互作用的模型，最后给出了应用该模型模拟液滴在固壁上铺展与滑行的结果。

2. 理论与模型

2.1. 格子Boltzmann方法

Boltzmann方程描述了由于其分子的流动和碰撞而导致的流体状态的演化。LBM是通过对Boltzmann方程的相空间进行离散化得到的，其中大量的分子速度被一组有限的速度矢量(晶格)代替 [1] 。LBM模型包含三要素：格子结构、流体粒子的离散速度集合、演化方程，其描述了流体粒子分布函数在固定格子上的演化过程

$f_i\left(x+e_i{\delta }_i,t+{\delta }_i\right)-f_i\left(x,t\right)={\Omega }_i\left(x,t\right)$ (1)

式中 $x$ 是格子上的格点， $\left\{e_i:i=1,2,\cdots ,m\right\}$ 是流体粒子的离散速度集合， ${\delta }_t$ 是时间增量，t是当前时间步， $f_i$ 是以速度 $e_i$ 运动的速度分布函数， $$ 是碰撞算子，表示流体粒子间的碰撞对速度分布函数的影响。流体的密度和宏观速度由离散分布函数的和得到

$\rho =m{\displaystyle {\sum }_i{f}_i}$ , $\rho u=m{\displaystyle {\sum }_i{e}_i}{f}_i$ (2)

m为流体粒子质量。

2.2. 多相或多组分格子Boltzmann模型

常用的多相或多组分格子Boltzmann模型为颜色模型、自由能模型和Shan-Chen模型。由于Shan-Chen模型具有能模拟大密度差多相流的优点，本文采用Shan-Chen模型。Shan-Chen模型考虑了具有不同分子质量的任意数量组分的密度分布函数。粒子之间的相互作用通过一组势包含在动力学中。第k组分的LB方程可写成以下形式：

$f_i^k\left(x+e_i{\delta }_i,t+{\delta }_i\right)-f_i^k\left(x,t\right)=-\frac{f_i^k\left(x,t\right)-f_i^{k\left(eq\right)}\left(x,t\right)}{{\tau }_k}$ (3)

式中， $f_i^k\left(x,t\right)$ 是第i个速度方向上第k组分在位置 $x$ 和时间t处的数密度分布， ${\delta }_t$ 是时间增量。在右手边， ${\tau }_k$ 是组分k的弛豫时间，以晶格单位表示， $f_i^{k\left(eq\right)}\left(x,t\right)$ 是相应的平衡分布函数。对于二维9-速度LB模型(D2Q9)，(其中D是维数，Q是速度方向数)， $f_i^{k\left(eq\right)}\left(x,t\right)$ 具有以下形式

$f_0^{k\left(eq\right)}={\alpha }_k{n}_k-\frac{2}{3}{n}_k{u}_k^{eq}\cdot u_k^{eq}$ (4)

$f_i^{k\left(eq\right)}=\frac{\left(1-{\alpha }_k\right)n_k}{5}+\frac{1}{3}{n}_k\left(e_i\cdot u_k^{eq}\right)+\frac{1}{2}{n}_k{\left(e_i\cdot u_k^{eq}\right)}^2-\frac{1}{6}{n}_k{u}_k^{eq}{u}_k^{eq}i=1,2,3,4$ (5)

$f_i^{k\left(eq\right)}=\frac{\left(1-{\alpha }_k\right)n_k}{20}+\frac{1}{12}{n}_k\left(e_i\cdot u_k^{eq}\right)+\frac{1}{8}{n}_k{\left(e_i\cdot u_k^{eq}\right)}^2-\frac{1}{24}{n}_k{u}_k^{eq}{u}_k^{eq}i=5,6,7,8$ (6)

在上述方程中，离散速度 $e_i$ 取为

$e_i=\{\begin{array}{l}0i=0\\ \left(\cos \frac{\left(i-1\right)\text{π}}{2},\sin \frac{\left(i-1\right)\text{π}}{2}\right)i=1,2,3,4\\ \sqrt{\text{2}}\left(\text{cos}\left[\frac{\left(i-5\right)\text{π}}{2}+\frac{\text{π}}{4}\right],\sin \left[\frac{\left(i-5\right)\text{π}}{2}+\frac{\text{π}}{4}\right]\right)i=5,6,7,8\end{array}$ (7)

${\alpha }_k$ 是一个自由参数，它与纯组分k区域的声速有关，其关系式为 ${\left(c_s^k\right)}^2=\frac{3}{5}\left(1-{\alpha }_k\right)$ ； $n_k={\displaystyle {\sum }_i{f}_i^k}$ 是k组分的数密度。k组分的质量密度 ${\rho }_k$ 定义为 ${\rho }_k=m_k{n}_k=m_k{\displaystyle {\sum }_i{f}_i^k}$ ，组分k的流体速度 $u_k$ 定义为 ${\rho }_k{u}_k=m_k{\displaystyle {\sum }_i{e}_i}{f}_i^k$ ，其中 $m_k$ 是k组分的分子质量。参数 $u_k^{eq}$ 由如下关系式确定：

${\rho }_k{u}_k^{eq}={\rho }_k{u}^{\prime }+{\tau }_k{F}_k$ (8)

其中， $F_k=F_{1k}+F_{2k}+F_{3k}$ 是作用于第k组分的总力，包括流体–流体相互作用 $F_{1k}$ 、外力 $F_{2k}$ 、流体–固体相互作用 $F_{3k}$ 。为使每个碰撞的动量守恒， $u^{\prime }$ 必须满足关系

$u^{\prime }=\left({\displaystyle \underset{k=1}{\overset{s}{\sum }}\frac{{\rho }_k{u}_k}{{\tau }_k}}\right)/\left({\displaystyle \underset{k=1}{\overset{s}{\sum }}\frac{{\rho }_k}{{\tau }_k}}\right)$ (9)

假设第k组分在 $x$ 点和第 $\bar{k}$ 组分在 $x^{\prime }$ 点之间的粒子相互作用力与它们的“有效数密度”的乘积 ${\Psi }_k\left(n_k\right)$ 成正比， ${\Psi }_k\left(n_k\right)$ 定义为局部数密度的函数。在k组分点 $x$ 处的总流体-流体相互作用力为

$F_{1k}\left(x\right)=-{\Psi }_k\left(x\right){\displaystyle \underset{x^{\prime }}{\sum }{\displaystyle \underset{k=1}{\overset{s}{\sum }}{G}_{k\bar{k}}}}\left(x,x^{\prime }\right){\Psi }_{\bar{k}}\left(x^{\prime }\right)\left(x-x^{\prime }\right)$ (10)

其中 $G_{k\bar{k}}\left(x,x^{\prime }\right)$ 满足 $G_{k\bar{k}}\left(x,x^{\prime }\right)=G_{k\bar{k}}\left(x^{\prime },x\right)$ ，而 ${\Psi }_k\left(x\right)$ 是 $x$ 处数密度 $n_k$ 的一个函数。

$G_{k\bar{k}}\left(x,x^{\prime }\right)=\{\begin{array}{l}{g}_{k\bar{k}},\left|x-x^{\prime }\right|=1\\ g_{k\bar{k}}/4,\left|x-x^{\prime }\right|=\sqrt{2}\\ 0,其它\end{array}$ (11)

这里， $g_{k\bar{k}}$ 是组分k和 $\bar{k}$ 粒子之间作用势的强度。本研究中，有效数密度 ${\Psi }_k\left(n_k\right)$ 取为 $n_k$ 。

常数体积力的作用可以简单地描述为

$F_{2k}={\rho }_{k}g=m_k{n}_{k}g$ (12)

式中 $g$ 是单位质量的体积力。

2.3. 流体与固壁之间的相互作用

在流体/固体界面，固壁被视为具有恒定数密度的相。流体与壁面之间的作用力描述为

$F_{3k}\left(x\right)=-n_k\left(x\right){\displaystyle \underset{x^{\prime }}{\sum }{g}_{kw}}{n}_w\left(x^{\prime }\right)\left(x^{\prime }-x\right)$ (13)

式中， $n_w$ 是墙的数量密度，在墙处为常数，在其他地方为零， $g_{kw}$ 是组分k与墙之间的相互作用强度。 $g_{kw}$ 对非润湿流体为正，对润湿流体为负。通过调节，我们可以得到不同的润湿性。

实施Chapman-Enskog展开步骤可以获得流体混合物作为单一流体的的连续性和动量方程：

$\frac{\partial \rho }{\partial t}+\nabla \cdot \left(\rho u\right)=0$ (14)

$\rho \left[\frac{\partial u}{\partial t}+(u\cdot \nabla )u\right]=-\nabla p+\nabla \cdot \left[\rho \nu \left(\nabla u+u\nabla \right)\right]+\rho g$ (15)

式中， $\rho ={\displaystyle {\sum }_k{\rho }_k}$ 是流体混合物的总密度，总流体速度 $u$ 由 $\rho u={\displaystyle {\sum }_k{\rho }_k{u}_k}+\frac{1}{2}{\displaystyle {\sum }_k{F}_k}$ 定义，请注意，引入流体-固体相互作用对宏观方程没有影响，因为 $F_{3k}$ 只存在于固体/流体界面。

3. 数值模拟

应用上述模型，模拟重力作用下液滴在水平固壁上铺展演化过程以及液滴在垂直壁面上的滑行过程。液滴铺展过程如图1所示，在初始时刻，液滴为半圆形，在重力作用下，液滴沿着固壁铺展，液滴高度逐渐减小，液滴与壁面的接触线逐渐增长。液滴滑行过程如图2所示，恒定的重力沿壁面方向作用于两相流体，密度大的液滴沿着壁面下滑，液滴逐渐变形。随着时间的推进，液滴变化逐渐减弱并趋于稳定。这些模拟结果与VOF模型的模拟结果一致 [3] 。

Figure 1. Evolution of droplet spreading on solid wall

图1. 液滴在壁面上的铺展演化过程

Figure 2. Evolution of droplet slipping on solid wall

图2. 液滴在壁面上的滑移演化过程

4. 总结

本文基于Shan-Chen模型建立了固壁上液滴动力学模型，并成功地实现了该算法，模拟了重力作用下液滴在固体表面上的铺展以及滑移动力学行为。

基金项目

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

参考文献

NOTES

^*通讯作者。