LBM笔记：Exact Difference Method 构建的源项

本文同步发表于知乎和个人主页.

我在阅读其他文献时，看到有文献提到 Kupershtokh 的一篇关于 LBM 外力项的文章¹，他提出的格式在其他文章中又被称为 Exact difference method（EDM）。

本想搜索这篇文章的PDF进行阅读，但最后只找到一个网站介绍了部分信息。这份粗略笔记的内容，都是基于其他文献^{2 3 4}中对该方法的介绍而整理的。

介绍

Boltzmann-BGK方程的近似

密度分布函数 $f$ 的 Boltzmann-BGK方程写作：

$$\frac{\partial f}{\partial t}+\xi \cdot \nabla f+\mathbf{a}\cdot {\nabla }_{\xi }f=-\frac{1}{{\tau }_f}\left[f-f^{(eq)}\right]\text{\hspace{1em}(1)}$$

其中：

$$f^{(eq)}=\frac{\rho }{(2\pi R_{g}T{)}^{D/2}}\exp \left(-\frac{(\xi -\mathbf{u}{)}^2}{2R_{g}T}\right)\text{\hspace{1em}(2)}$$

符号含义分别为

• $\xi$ : 粒子速度；
• $\rho$ : 流体密度；
• $\mathbf{u}$ : 流场宏观速度；
• $R_g$ : 气体常数；
• $T$ : 温度；
• $\mathbf{a}=\mathrm{d}\mathbf{u}/\mathrm{d}t$ : 流场加速度。

当Knudsen数较小时，Boltzmann方程的体力项可以基于$f^{(eq)}$进行近似，即
$${\nabla }_{\xi }f\approx {\nabla }_{\xi }{f}^{(eq)}$$

根据式(2)， $f^{(eq)}$ 是分子随机速度 $\mathbf{C}=\xi -\mathbf{u}$ 的指数函数，则：
$$\frac{\partial f^{(eq)}}{\partial \mathbf{u}}=-\frac{\partial f^{(eq)}}{\partial \xi }=-{\nabla }_{\xi }{f}^{(eq)}$$

若假设密度 $\rho$ 为常数，则 $\frac{\mathrm{d}\rho }{\mathrm{d}t}=0$，那么 $f^{(eq)}$ 的全导数为：

$$\begin{array}{rl}\frac{\mathrm{d}{f}^{(eq)}}{\mathrm{d}t}& =\frac{\partial f^{(eq)}}{\partial \mathbf{u}}\cdot \frac{\mathrm{d}\mathbf{u}}{\mathrm{d}t}+\frac{\partial f^{(eq)}}{\partial \rho }\cdot \frac{\mathrm{d}\rho }{\mathrm{d}t}\\ & \approx \frac{\partial f^{(eq)}}{\partial \mathbf{u}}\cdot \frac{\mathrm{d}\mathbf{u}}{\mathrm{d}t}=-\mathbf{a}\cdot {\nabla }_{\xi }{f}^{(eq)}\end{array}$$

也就是说，Kupershtokh 等¹²近似的 Boltzmann-BGK方程写作：

$$\frac{\partial f}{\partial t}+\xi \cdot \nabla f=-\frac{1}{{\tau }_f}\left[f-f^{(eq)}\right]+\frac{\mathrm{d}{f}^{(eq)}}{\mathrm{d}t}\text{\hspace{1em}(3)}$$

近似方程的离散化

将式(3)在$[t,t+{\delta }_t]$内沿着${\mathbf{e}}_{\alpha }$方向积分，可得

$$\begin{array}{rl}{f}_{\alpha }(\mathbf{x}+{\mathbf{e}}_{\alpha }{\delta }_t,t+{\delta }_t)-f_{\alpha }(\mathbf{x},t)=& -\frac{1}{\tau }\left[f_{\alpha }(\mathbf{x},t)-f_{\alpha }^{(eq)}(\rho ,{\mathbf{u}}^{\ast })\right]\\ & +\left[f_{\alpha }^{(eq)}({\mathbf{u}}^{\ast }(t+{\delta }_t))-f_{\alpha }^{(eq)}({\mathbf{u}}^{\ast }(t))\right]\end{array}\text{\hspace{1em}(4)}$$

其中：

• ${\delta }_t$ 为时间步长；
• $\tau ={\tau }_f/{\delta }_t$ 为无量纲松弛时间；
• $\mathbf{x}$ 为网格点位置，$\rho ,{\mathbf{u}}^{\ast }$ 分别为该点在 $t$ 时刻的密度和速度（与宏观速度 $\mathbf{u}$ 有区别）；
• $f_{\alpha },f_{\alpha }^{(eq)}$ 分别为 ${\mathbf{e}}_{\alpha }$ 方向的分布函数和平衡态；
  $$\begin{array}{rl}{f}_{\alpha }^{(eq)}& =w_{\alpha }\rho \cdot \left[1+\frac{{\mathbf{e}}_{\alpha }\cdot \mathbf{u}}{c_s^2}+\frac{({\mathbf{e}}_{\alpha }\cdot \mathbf{u}{)}^2}{2c_s^4}-\frac{|\mathbf{u}{|}^2}{2c_s^2}\right]\\ & =w_{\alpha }\rho \cdot \left[1+\frac{{\mathbf{e}}_{\alpha }\cdot \mathbf{u}}{c_s^2}+\frac{1}{2c_s^4}(\mathbf{uu}):({\mathbf{e}}_{\alpha }{\mathbf{e}}_{\alpha }-c_s^2[\mathbf{I}])\right]\end{array}\text{\hspace{1em}(4-1)}$$
  这里 $w_{\alpha }$ 为 ${\mathbf{e}}_{\alpha }$ 的权重系数， $[\mathbf{I}]$ 为单位矩阵；两个同样为m*n大小的矩阵 $[\mathbf{A}]$ 和 $[\mathbf{B}]$ ，则 $[\mathbf{A}]:[\mathbf{B}]={\sum }_{i=1}^m{\sum }_{j=1}^n{A}_{ij}{B}_{ij}$ ；
• $F_{\alpha }$ 为外力项，其表达式基于外力 $\mathbf{F}=\rho \frac{\mathrm{d}{\mathbf{u}}^{\ast }}{\mathrm{d}t}=\rho \mathbf{a}$ 进行构建。

$$F_{\alpha }=f_{\alpha }^{(eq)}({\mathbf{u}}^{\ast }(t+{\delta }_t))-f_{\alpha }^{(eq)}({\mathbf{u}}^{\ast }(t))$$

并且

$$\begin{array}{rl}{\mathbf{u}}^{\ast }(t)& =\left(\sum _{\alpha }{f}_{\alpha }{\mathbf{e}}_{\alpha }\right)/\left(\sum _{\alpha }{f}_{\alpha }\right)\\ {\mathbf{u}}^{\ast }(t+{\delta }_t)& ={\mathbf{u}}^{\ast }(t)+{\int }_t^{t+{\delta }_t}\frac{\mathrm{d}{\mathbf{u}}^{\ast }}{\mathrm{d}t}\mathrm{d}t\\ & ={\mathbf{u}}^{\ast }(t)+{\int }_t^{t+{\delta }_t}\frac{\mathbf{F}}{\rho }\mathrm{d}t\approx {\mathbf{u}}^{\ast }(t)+\frac{\mathbf{F}}{\rho }{\delta }_t\end{array}\text{\hspace{1em}(5)}$$

这意味着加速度 $\mathbf{a}=\mathbf{F}/\rho$ 在 $[t,t+{\delta }_t]$ 时间段内被视为常数。

Kupershtokh等在这篇文章²说明，在使用该源项时，流体宏观量的计算应表示为：

$$\rho =\sum _{\alpha }{f}_{\alpha },\rho \mathbf{u}=\sum _i{\mathbf{e}}_{\alpha }{f}_{\alpha }+\frac{{\delta }_t\mathbf{F}}{2}\text{\hspace{1em}(6)}$$

Li 等⁴通过 Chapman-Enskog 展开证明：该格式向宏观方程引入的误差项为：
$${\mathbf{Error}}_{\mathrm{EDM}}=-{\delta }_t^2\nabla \cdot \left(\frac{\mathbf{FF}}{4\rho }\right)\text{\hspace{1em}(7)}$$

EDM 源项的展开式

为便于分析，这里采用如下书写形式的平衡态：
$$f_{\alpha }^{(eq)}=w_{\alpha }\rho \cdot \left[1+\frac{{\mathbf{e}}_{\alpha }\cdot \mathbf{u}}{c_s^2}+\frac{1}{2c_s^4}(\mathbf{uu}):({\mathbf{e}}_{\alpha }{\mathbf{e}}_{\alpha }-c_s^2[\mathbf{I}])\right]$$

将其带入 EDM 源项，记 $\Delta \mathbf{u}=\frac{\mathbf{F}}{\rho }{\delta }_t$，整理得：

$$\begin{gathered} \begin{array}{rl}{F}_{\alpha }& =f_{\alpha }^{(eq)}({\mathbf{u}}^{\ast }(t+{\delta }_t))-f_{\alpha }^{(eq)}({\mathbf{u}}^{\ast }(t))\\ & =\rho w_{\alpha }\{\frac{{\mathbf{e}}_{\alpha }\cdot \Delta \mathbf{u}}{c_s^2}+ \\ \frac{1}{2c_s^4}[({\mathbf{u}}^{\ast }+\Delta \mathbf{u})({\mathbf{u}}^{\ast }+\Delta \mathbf{u})\\ & -{\mathbf{u}}^{\ast }{\mathbf{u}}^{\ast }]:({\mathbf{e}}_{\alpha }{\mathbf{e}}_{\alpha }-c_s^2[\mathbf{I}])\}\end{array} \end{gathered}$$

这里涉及矩阵 $({\mathbf{u}}^{\ast }+\Delta \mathbf{u})({\mathbf{u}}^{\ast }+\Delta \mathbf{u})$ 和 ${\mathbf{u}}^{\ast }{\mathbf{u}}^{\ast }$ 的相减，由于：
$$({\mathbf{u}}^{\ast }+\Delta \mathbf{u})({\mathbf{u}}^{\ast }+\Delta \mathbf{u})=(u_i^{\ast }+\Delta u_i)(u_j^{\ast }+\Delta u_j),{\mathbf{u}}^{\ast }{\mathbf{u}}^{\ast }=u_i^{\ast }{u}_j^{\ast }$$

其中 $u_i$ 和 $\Delta u_i=F_i{\delta }_t/\rho$ 分别为 $\mathbf{u}$ 和 $\Delta \mathbf{u}$ 在 i 轴上的分量，所以：
$$\begin{array}{rl}& ({\mathbf{u}}^{\ast }+\Delta \mathbf{u})({\mathbf{u}}^{\ast }+\Delta \mathbf{u})-{\mathbf{u}}^{\ast }{\mathbf{u}}^{\ast }\\ =& (u_i^{\ast }+\Delta u_i)(u_j^{\ast }+\Delta u_j)-u_i^{\ast }{u}_j^{\ast }\\ =& \frac{{\delta }_t}{\rho }(u_i^{\ast }{F}_j+u_j^{\ast }{F}_i)+\frac{{{\delta }_t}^2}{{\rho }^2}{F}_i{F}_j\\ =& \frac{{\delta }_t}{\rho }(u_i^{\ast }+\frac{F_i{\delta }_t}{2\rho })F_j+\frac{{\delta }_t}{\rho }(u_j^{\ast }+\frac{F_j{\delta }_t}{2\rho })F_i\end{array}$$

这里令 ${\mathbf{v}}_{\mathrm{EDM}}={\mathbf{u}}^{\ast }+\frac{\mathbf{F}{\delta }_t}{2\rho }$ ，则 EDM 源项的展开式化简为：

$$\begin{array}{rl}{F}_{\alpha }& ={\delta }_t{w}_{\alpha }\left[\frac{{\mathbf{e}}_{\alpha }\cdot \mathbf{F}}{c_s^2}+\frac{1}{2c_s^4}\left({\mathbf{v}}_{\mathrm{EDM}}\mathbf{F}+\mathbf{F}{\mathbf{v}}_{\mathrm{EDM}}\right):({\mathbf{e}}_{\alpha }{\mathbf{e}}_{\alpha }-c_s^2[\mathbf{I}])\right]\end{array}\text{\hspace{1em}(8)}$$
上式即为EDM源项的展开结果。

---

1. Kupershtokh, A. L. New method of incorporating a body force term into the lattice Boltzmann equation. in Proceeding of the 5th international EHD workshop 241–246 (Poitiers, France, 2004). ↩︎ ↩︎

2. A.L. Kupershtokh, D.A. Medvedev, & D.I. Karpov (2009). On equations of state in a lattice Boltzmann method. Computers & Mathematics with Applications, 58(5), 965-974. DOI:10.1016/j.camwa.2009.02.024. ↩︎ ↩︎ ↩︎

3. Q. Li, K.H. Luo, Q.J. Kang, Y.L. He, Q. Chen, & Q. Liu (2016). Lattice Boltzmann methods for multiphase flow and phase-change heat transfer. Progress in Energy and Combustion Science, 52, 62-105. DOI:10.1016/j.pecs.2015.10.001. ↩︎

4. Li, Q., Zhou, P., & Yan, H. (2016). Revised Chapman-Enskog analysis for a class of forcing schemes in the lattice Boltzmann method. Physical Review E, 94, 043313. DOI:10.1103/PhysRevE.94.043313. ↩︎ ↩︎