Literature Review : 简化的格子 Boltzmann 方法

注：该文章同步发在知乎专栏

这里所说的 “简化的格子 Boltzmann 方法” 其实是 Simplified and Highly Stable Lattice Boltzmann Method (SHSLBM)(¹,²,³,⁴)。该方法的特色是仅使用平衡态分布函数、宏观密度、宏观速度进行 LBM 计算，并保证无条件稳定。

下面的内容只是我对相关文献进行简单阅读后的部分整理

单松弛格子 Boltzmann 方法

基本方程

针对流场的 LBGK 演化方程为

$$f_{\alpha }(\mathbf{r}+{\mathbf{c}}_{\alpha }+\Delta t,t+\Delta t)-f_{\alpha }(\mathbf{r},t)=-\frac{1}{\tau }\left[f_{\alpha }(\mathbf{r},t)-f_{\alpha }^{(eq)}\right]\text{\hspace{1em}(1)}$$

其中 $f_{\alpha }(\mathbf{r},t)$ 为 $t$ 时刻在 $\mathbf{r}$ 处沿着 ${\mathbf{c}}_{\alpha }$ 方向的密度分布函数， $\tau$ 为 BGK 碰撞项的松弛时间。 $f_{\alpha }^{(eq)}$ 为 ${\mathbf{c}}_{\alpha }$ 方向的平衡态分布

$$f_{\alpha }^{(eq)}=\rho w_{\alpha }\left[1+\frac{{\mathbf{c}}_{\alpha }\cdot \mathbf{u}}{c_s^2}+\frac{({\mathbf{c}}_{\alpha }\cdot \mathbf{u}{)}^2}{2c_s^4}-\frac{\mathbf{u}\cdot \mathbf{u}}{2c_s^2}\right]\text{\hspace{1em}(2)}$$

式 (2) 中 $\rho$ 为流体密度， $\mathbf{u}=(u_x,u_y,u_z{)}^T$ 为流体宏观速度， $w_{\alpha }$ 为 ${\mathbf{c}}_{\alpha }$ 的权重系数， $c_s$ 为格子声速。 $w_{\alpha }$ 、 ${\mathbf{c}}_{\alpha }$ 和 $c_s$ 的数值由选取的离散速度模型决定。

对于 LBM 而言，宏观量表示为

$$\rho =\sum _{\alpha }{f}_{\alpha }=\sum _{\alpha }{f}_{\alpha }^{(eq)},\rho \mathbf{u}=\sum _{\alpha }{\mathbf{c}}_{\alpha }{f}_{\alpha }=\sum _{\alpha }{\mathbf{c}}_{\alpha }{f}_{\alpha }^{(eq)}\text{\hspace{1em}(3)}$$

Chapman-Enskog 展开

Boltzmann 方程和 LBGK 方程都可以通过多尺度展开还原为 Navier-Stokes 方程。这里以 Chapman-Enskog 为例：
$$f_{\alpha }=f_{\alpha }^{(eq)}+ϵf_{\alpha }^{(1)}+{ϵ}^2{f}_{\alpha }^{(2)}+...,\frac{\partial }{\partial t}=ϵ\frac{\partial }{\partial t_0}+{ϵ}^2\frac{\partial }{\partial t_1},\frac{\partial }{\partial \mathbf{r}}=ϵ\frac{\partial }{\partial {\mathbf{r}}_0}$$

上述展开保留至 ${ϵ}^1$ ，即可还原为⁵
$$\frac{\partial \rho }{\partial t}+\nabla \cdot \left(\sum _{\alpha }{\mathbf{c}}_{\alpha }{f}_{\alpha }^{(eq)}\right)=0\text{\hspace{1em}(4)}$$

和
$$\frac{\partial \rho \mathbf{u}}{\partial t}+\nabla \cdot \Pi =0\text{\hspace{1em}(5)}$$

式 (5) 中张量 $\Pi$ 表示为
$${\Pi }_{\beta \gamma }=\sum _{\alpha }({\mathbf{c}}_{\alpha }{)}_{\beta }({\mathbf{c}}_{\alpha }{)}_{\gamma }\left[f_{\alpha }^{(eq)}+\left(1-\frac{1}{2\tau }\right)ϵf_{\alpha }^{(1)}\right]$$

由于只保留到 ${ϵ}^1$ ，所以 $ϵf_{\alpha }^{(1)}\approx f_{\alpha }-f_{\alpha }^{(eq)}=f_{\alpha }^{(neq)}$ 近似为非平衡态分布函数。此时
$${\Pi }_{\beta \gamma }=\sum _{\alpha }({\mathbf{c}}_{\alpha }{)}_{\beta }({\mathbf{c}}_{\alpha }{)}_{\gamma }\left[f_{\alpha }^{(eq)}+\left(1-\frac{1}{2\tau }\right)f_{\alpha }^{(neq)}\right]\text{\hspace{1em}(6)}$$

并且存在如下关系
$$f_{\alpha }^{(neq)}\approx ϵf_{\alpha }^{(1)}=-\tau \left[\frac{\partial }{\partial t}+{\mathbf{c}}_{\alpha }\cdot \nabla \right]f_{\alpha }^{(eq)}\Delta t\text{\hspace{1em}(7)}$$

令 $\mathrm{D}=\frac{\partial }{\partial t}+{\mathbf{c}}_{\alpha }\cdot \nabla$ ，则 $f_{\alpha }^{(neq)}\approx -\tau \cdot (\Delta t)\cdot \mathrm{D}{f}_{\alpha }^{(eq)}$。

流体运动粘度 $\nu$ 和松弛时间 $\tau$ 的关系为
$$\nu =\left(\tau -\frac{1}{2}\right)c_s^2\Delta t$$

SHSLBM

式 (7) 说明：在保留至 ${ϵ}^1$ 的情况下， $f_{\alpha }^{(neq)}$ 和 $f_{\alpha }^{(eq)}$ 存在联系。考虑到式 (1) 中由 $\frac{-1}{\tau }$ 松弛的部分即为 $f_{\alpha }^{(neq)}$ ，因此式 (1) 理应可以被写作完全由 $f_{\alpha }^{(eq)}$ 进行表示的形式。这便是 SHSLBM 的核心思想。

基本算法

SHSLBM 的单个时间步迭代分为预报和校正两部分，具体写作：

1. 预报步：计算 $(\mathbf{r},t)$ 的预估密度 ${\rho }^{\ast }$ 和预估速度 ${\mathbf{u}}^{\ast }$

$${\rho }^{\ast }=\sum _{\alpha }{f}_{\alpha }^{(eq)}(\mathbf{r}-{\mathbf{c}}_{\alpha }\Delta t,t-\Delta t),{\rho }^{\ast }{\mathbf{u}}^{\ast }=\sum _{\alpha }{\mathbf{c}}_{\alpha }{f}_{\alpha }^{(eq)}(\mathbf{r}-{\mathbf{c}}_{\alpha }\Delta t,t-\Delta t)\text{\hspace{1em}(8)}$$

2. 校正步：根据 ${\rho }^{\ast }$ 和 ${\mathbf{u}}^{\ast }$ 得出校正值

$$\{\begin{array}{l}\rho ={\rho }^{\ast }\\ \rho \mathbf{u}={\rho }^{\ast }{\mathbf{u}}^{\ast }-\left(1-\frac{1}{\tau }\right)\sum _{\alpha }{\mathbf{c}}_{\alpha }[f_{\alpha }^{(neq)}(\mathbf{r}+\frac{\Delta t}{2}{\mathbf{c}}_{\alpha },t-\frac{\Delta t}{2})-f_{\alpha }^{(neq)}(\mathbf{r}-\frac{\Delta t}{2}{\mathbf{c}}_{\alpha },t-\frac{\Delta t}{2})]\end{array}\text{\hspace{1em}(9-1)}$$

式 (9-1) 中非平衡态的计算通过对式 (7) 的差分实现，即：

$$f_{\alpha }^{(neq)}(\mathbf{r}-\frac{\Delta t}{2}{\mathbf{c}}_{\alpha },t-\frac{\Delta t}{2})=-\tau \left[f_{\alpha }^{(eq)}(\mathbf{r},t)-f_{\alpha }^{(eq)}(\mathbf{r}-{\mathbf{c}}_{\alpha }\Delta t,t-\Delta t)\right]+O((\Delta t{)}^3)$$

同理：

$$f_{\alpha }^{(neq)}(\mathbf{r}+\frac{\Delta t}{2}{\mathbf{c}}_{\alpha },t-\frac{\Delta t}{2})=-\tau \left[f_{\alpha }^{(eq)}(\mathbf{r}+{\mathbf{c}}_{\alpha }\Delta t,t)-f_{\alpha }^{(eq)}(\mathbf{r},t-\Delta t)\right]+O((\Delta t{)}^3)$$

在这两条式子中， $f_{\alpha }^{(eq)}(\mathbf{r},t)$ 和 $f_{\alpha }^{(eq)}(\mathbf{r}+{\mathbf{c}}_{\alpha }\Delta t,t)$ 都是未知的。参照Shu等⁶的做法，对上面两条公式中的这两个量使用预报步数据进行近似，即

$$f_{\alpha }^{(eq)}(\mathbf{r},t)=f_{\alpha }^{(eq)}({\rho }^{\ast }(\mathbf{r},t),{\mathbf{u}}^{\ast }(\mathbf{r},t))$$

式 (9-1) 中 $\rho \mathbf{u}$ 的计算还可以进一步化简。注意到，对于 $f_{\alpha }^{(neq)}(\mathbf{r}-\frac{\Delta t}{2}{\mathbf{c}}_{\alpha },t-\frac{\Delta t}{2})$ ，有

$$\begin{array}{rl}\sum _{\alpha }{\mathbf{c}}_{\alpha }{f}_{\alpha }^{(neq)}(\mathbf{r}-\frac{\Delta t}{2}{\mathbf{c}}_{\alpha },t-\frac{\Delta t}{2})& =-\tau [\sum _{\alpha }{\mathbf{c}}_{\alpha }{f}_{\alpha }^{(eq)}(\mathbf{r},t)\\ & -\sum _{\alpha }{\mathbf{c}}_{\alpha }{f}_{\alpha }^{(eq)}(\mathbf{r}-{\mathbf{c}}_{\alpha }\Delta t,t-\Delta t)]+O((\Delta t{)}^3)\\ & =-\tau ({\rho }^{\ast }{\mathbf{u}}^{\ast }-{\rho }^{\ast }{\mathbf{u}}^{\ast })+O((\Delta t{)}^3)\\ & =O((\Delta t{)}^3)\end{array}$$

所以式 (9-1) 可进一步化简为式 (9-2)

$$\{\begin{array}{l}\rho ={\rho }^{\ast }\\ \rho \mathbf{u}={\rho }^{\ast }{\mathbf{u}}^{\ast }-\left(1-\frac{1}{\tau }\right)\sum _{\alpha }{\mathbf{c}}_{\alpha }{f}_{\alpha }^{(neq)}(\mathbf{r}+\frac{\Delta t}{2}{\mathbf{c}}_{\alpha },t-\frac{\Delta t}{2})\end{array}\text{\hspace{1em}(9-2)}$$

对应的宏观流程

事实上，式 (8) 和式 (9-1), (9-2) 的预报和校正存在着对应的宏观偏微分方程。下面简要说明一下。

预报步

式 (8) 对应的宏观方程为

$$\{\begin{array}{l}\frac{\partial \rho }{\partial t}+\nabla \cdot \left(\sum _{\alpha }{\mathbf{c}}_{\alpha }{f}_{\alpha }^{(eq)}\right)=0\\ \frac{\partial \rho \mathbf{u}}{\partial t}+\nabla \cdot \left[\sum _{\alpha }({\mathbf{c}}_{\alpha }{)}_{\beta }({\mathbf{c}}_{\alpha }{)}_{\gamma }{f}_{\alpha }^{(eq)}+\frac{1}{2\tau }\sum _{\alpha }({\mathbf{c}}_{\alpha }{)}_{\beta }({\mathbf{c}}_{\alpha }{)}_{\gamma }{f}_{\alpha }^{(neq)}\right]=0\end{array}\text{\hspace{1em}(10)}$$

首先，对 $f_{\alpha }^{(eq)}(\mathbf{r}-{\mathbf{c}}_{\alpha }\Delta t,t-\Delta t)$ 在 $(\mathbf{r},t)$ 进行展开，并代入式 (7) ，得到：

$$\begin{array}{rl}{f}_{\alpha }^{(eq)}(\mathbf{r}-{\mathbf{c}}_{\alpha }\Delta t,t-\Delta t)& =f_{\alpha }^{(eq)}(\mathbf{r},t)-(\Delta t)\mathrm{D}{f}_{\alpha }^{(eq)}(\mathbf{r},t)\\ & +\frac{(\Delta t{)}^2}{2}{\mathrm{D}}^2{f}_{\alpha }^{(eq)}(\mathbf{r},t)+O((\Delta t{)}^3)\\ & =f_{\alpha }^{(eq)}(\mathbf{r},t)-(\Delta t)\mathrm{D}{f}_{\alpha }^{(eq)}(\mathbf{r},t)\\ & -\frac{\Delta t}{2\tau }\mathrm{D}{f}_{\alpha }^{(neq)}(\mathbf{r},t)+O((\Delta t{)}^3)\end{array}\text{\hspace{1em}(11)}$$

将式 (11) 代入式 (8) 的第一条公式，有：
$$\begin{array}{rl}{\rho }^{\ast }& =\sum _{\alpha }{f}_{\alpha }^{(eq)}(\mathbf{r},t)-(\Delta t)\sum _{\alpha }\mathrm{D}{f}_{\alpha }^{(eq)}(\mathbf{r},t)\\ & -\frac{\Delta t}{2\tau }\sum _{\alpha }\mathrm{D}{f}_{\alpha }^{(neq)}(\mathbf{r},t)+O((\Delta t{)}^3)\end{array}\text{\hspace{1em}(12)}$$

式 (12) 右侧第一项与 ${\rho }^{\ast }$ 相等，右侧第三项根据非平衡态的守恒律可知其为 0 。所以在消去式 (12) 中的 $\Delta t$ 后， ${\rho }^{\ast }$ 的计算与式 (10) 中第一条公式对应，并具有 $O((\Delta t{)}^2)$ 精度。

将式 (11) 代入式 (8) 的第二条公式，整理得：
$$\begin{array}{rl}{\rho }^{\ast }{\mathbf{u}}^{\ast }& =\sum _{\alpha }{\mathbf{c}}_{\alpha }{f}_{\alpha }^{(eq)}(\mathbf{r},t)-(\Delta t)\sum _{\alpha }{\mathbf{c}}_{\alpha }\mathrm{D}{f}_{\alpha }^{(eq)}(\mathbf{r},t)\\ & -\frac{\Delta t}{2\tau }\sum _{\alpha }{\mathbf{c}}_{\alpha }\mathrm{D}{f}_{\alpha }^{(neq)}(\mathbf{r},t)+O((\Delta t{)}^3)\\ & =\sum _{\alpha }{\mathbf{c}}_{\alpha }{f}_{\alpha }^{(eq)}(\mathbf{r},t)-(\Delta t)[\frac{\partial }{\partial t}\sum _{\alpha }{\mathbf{c}}_{\alpha }{f}_{\alpha }^{(eq)}(\mathbf{r},t)+\\ & \nabla \cdot \sum _{\alpha }({\mathbf{c}}_{\alpha }{)}_{\beta }({\mathbf{c}}_{\alpha }{)}_{\gamma }\left(f_{\alpha }^{(eq)}(\mathbf{r},t)+\frac{1}{2\tau }{f}_{\alpha }^{(neq)}(\mathbf{r},t)\right)]\\ & -\frac{\Delta t}{2\tau }\sum _{\alpha }\frac{\partial }{\partial t}{\mathbf{c}}_{\alpha }{f}_{\alpha }^{(neq)}(\mathbf{r},t)+O((\Delta t{)}^3)\end{array}\text{\hspace{1em}(13)}$$

式 (13) 右侧第一项与 ${\rho }^{\ast }{\mathbf{u}}^{\ast }$ 相等，右侧 $\sum _{\alpha }\frac{\partial }{\partial t}{\mathbf{c}}_{\alpha }{f}_{\alpha }^{(neq)}(\mathbf{r},t)$ 项根据非平衡态的守恒律可知其为 0 。所以在消去式 (13) 中的 $\Delta t$ 后， ${\rho }^{\ast }{\mathbf{u}}^{\ast }$ 的计算与式 (10) 中第二条公式对应，并具有 $O((\Delta t{)}^2)$ 精度。

校正步

式 (9) 对应的宏观方程为
$$\{\begin{array}{l}\frac{\partial \rho }{\partial t}=0\\ \frac{\partial \rho \mathbf{u}}{\partial t}+\nabla \cdot \left[\left(1-\frac{1}{2\tau }\right)\sum _{\alpha }({\mathbf{c}}_{\alpha }{)}_{\beta }({\mathbf{c}}_{\alpha }{)}_{\gamma }{f}_{\alpha }^{(neq)}\right]=0\end{array}\text{\hspace{1em}(14)}$$

这里主要介绍一下动量方程的还原。对 $f_{\alpha }^{(neq)}(\mathbf{r}\pm \frac{\Delta t}{2}{\mathbf{c}}_{\alpha },t-\frac{\Delta t}{2})$ 进行 Taylor 展开，得到

$$\begin{array}{rl}{f}_{\alpha }^{(neq)}(\mathbf{r}\pm \frac{\Delta t}{2}{\mathbf{c}}_{\alpha },t-\frac{\Delta t}{2})& =f_{\alpha }^{(neq)}(\mathbf{r},t-\frac{\Delta t}{2})\\ & \pm \frac{\Delta t}{2}({\mathbf{c}}_{\alpha }\cdot \nabla )f_{\alpha }^{(neq)}(\mathbf{r},t-\frac{\Delta t}{2})\\ & +\frac{(\Delta t{)}^2}{8}({\mathbf{c}}_{\alpha }\cdot \nabla {)}^2{f}_{\alpha }^{(neq)}(\mathbf{r},t-\frac{\Delta t}{2})+O((\Delta t{)}^3)\end{array}$$

这里我们以未化简的式 (9-1) 为例。在式 (9-1) 的动量校正步中，有：

$$\begin{array}{rl}& \left(1-\frac{1}{\tau }\right)\frac{1}{\Delta t}\sum _{\alpha }{\mathbf{c}}_{\alpha }[f_{\alpha }^{(neq)}(\mathbf{r}+\frac{\Delta t}{2}{\mathbf{c}}_{\alpha },t-\frac{\Delta t}{2})-f_{\alpha }^{(neq)}(\mathbf{r}-\frac{\Delta t}{2}{\mathbf{c}}_{\alpha },t-\frac{\Delta t}{2})]\\ & =\left(1-\frac{1}{\tau }\right)\sum _{\alpha }{\mathbf{c}}_{\alpha }({\mathbf{c}}_{\alpha }\cdot \nabla )f_{\alpha }^{(neq)}(\mathbf{r},t-\frac{\Delta t}{2})+O((\Delta t{)}^2)\\ & =\nabla \cdot \left[\left(1-\frac{1}{\tau }\right)\sum _{\alpha }({\mathbf{c}}_{\alpha }{)}_{\beta }({\mathbf{c}}_{\alpha }{)}_{\gamma }{f}_{\alpha }^{(neq)}(\mathbf{r},t-\frac{\Delta t}{2})\right]+O((\Delta t{)}^2)\end{array}$$

此时联立式 (13) ，即可还原式 (14) 中的动量方程。

SHSLBM 的稳定性分析

这里令 $\mathbf{Y}=(y_1,y_2,y_3,y_4{)}^T$。其中

$$y_1=\rho ,y_2=\rho u_x,y_3=\rho u_y,y_4=\rho u_z$$

将这些符号代换到式 (2) 的 $f_{\alpha }^{(eq)}$ 中，得：
$$\begin{array}{rl}{f}_{\alpha }^{(eq)}& =y_{1\alpha }{w}_{\alpha }+\frac{w_{\alpha }}{c_s^2}(c_{\alpha x}{y}_{2\alpha }+c_{\alpha y}{y}_{3\alpha }+c_{\alpha z}{y}_{4\alpha })\\ & +\frac{w_{\alpha }}{2y_{1\alpha }{c}_s^4}(c_{\alpha x}{y}_{2\alpha }+c_{\alpha y}{y}_{3\alpha }+c_{\alpha z}{y}_{4\alpha }{)}^2\\ & -\frac{w_{\alpha }}{2y_{1\alpha }{c}_s^2}(y_{2\alpha }^2+y_{3\alpha }^2+y_{4\alpha }^2)\end{array}\text{\hspace{1em}(15)}$$
where the subscript $\alpha$ denotes the quantities at a space level of $(\mathbf{r}-{\mathbf{c}}_{\alpha }\Delta t)$.²

预报步

对于第 $n$ 步的结果 ${\mathbf{Y}}^n$，式 (8) 的预报步写作 ${\mathbf{Y}}^{\ast }=\mathbf{G}({\mathbf{Y}}^n)$。

为便于进行诺依曼稳定性分析，对上式进行求导，有
$$\mathrm{d}{\mathbf{Y}}^{\ast }={\left[\frac{\partial {\mathbf{Y}}^{\ast }}{\partial y_1}\right]}^n\mathrm{d}{y}_1^n+{\left[\frac{\partial {\mathbf{Y}}^{\ast }}{\partial y_2}\right]}^n\mathrm{d}{y}_2^n+{\left[\frac{\partial {\mathbf{Y}}^{\ast }}{\partial y_3}\right]}^n\mathrm{d}{y}_3^n+{\left[\frac{\partial {\mathbf{Y}}^{\ast }}{\partial y_4}\right]}^n\mathrm{d}{y}_4^n$$

上标 $n$ 代表第 $n$ 个时间步。由于 $\mathrm{d}{\mathbf{Y}}^{\ast }=(\mathrm{d}{y}_1^{\ast },\mathrm{d}{y}_2^{\ast },\mathrm{d}{y}_3^{\ast },\mathrm{d}{y}_4^{\ast }{)}^T$ ，所以，对于预报步结果的每个分量而言有： $\mathrm{d}{y}_{\lambda }^{\ast }={\sum }_{\kappa =1}^4{\left[\frac{\partial y_{\lambda }^{\ast }}{\partial y_{\kappa }}\right]}^n\mathrm{d}{y}_{\kappa }^n$ （$\lambda ,\kappa \in [1,2,3,4]$）。

对于本文列出的方程组， ${\left[\frac{\partial y_{\lambda }^{\ast }}{\partial y_{\kappa }}\right]}^n$ 多达 16 个。因此下文仅拿 ${\left[\frac{\partial y_1^{\ast }}{\partial y_1}\right]}^n$ 举例进行计算。

$$\begin{array}{rl}{\left[\frac{\partial y_1^{\ast }}{\partial y_1}\right]}^n& =\frac{\partial }{\partial y_1}\sum _{\alpha }{f}_{\alpha }^{(eq)}\\ & =\sum _{\alpha }[w_{\alpha }-\frac{w_{\alpha }}{2(y_{1\alpha }^n{)}^2{c}_s^4}(c_{\alpha x}{y}_{2\alpha }^n+c_{\alpha y}{y}_{3\alpha }^n+c_{\alpha z}{y}_{4\alpha }^n{)}^2\\ & +\frac{w_{\alpha }}{2(y_{1\alpha }^n{)}^2{c}_s^2}((y_{2\alpha }^n{)}^2+(y_{3\alpha }^n{)}^2+(y_{4\alpha }^n{)}^2)]\end{array}\text{\hspace{1em}(16)}$$

假设式 (16) 中的 $\mathrm{d}{y}_i^{\ast }$, $\mathrm{d}{y}_i^n$, $y_{i\alpha }^n$ 均可被展开为如下形式（$i=[1,2,3,4]$，$j=\sqrt{-1}$）

$$\mathrm{d}{y}_i^{\ast }={\overline{\mathrm{d}{y}_i}}^{\ast }\exp (j\mathbf{k}\cdot \mathbf{r}),\mathrm{d}{y}_i^n={\overline{\mathrm{d}{y}_i}}^n\exp (j\mathbf{k}\cdot \mathbf{r}),y_{i\alpha }^n=y_i^n\exp (j\mathbf{k}\cdot \mathbf{r})$$

代入到式 (16) 和 $\mathrm{d}{\mathbf{Y}}^{\ast }$ 表达式，得

$$\begin{array}{rl}{\left[\frac{\partial y_1^{\ast }}{\partial y_1}\right]}^n& =\frac{\partial }{\partial y_1}\sum _{\alpha }{f}_{\alpha }^{(eq)}\\ & =\sum _{\alpha }[w_{\alpha }-\frac{w_{\alpha }}{2(y_1^n{)}^2{c}_s^4}(c_{\alpha x}{y}_2^n+c_{\alpha y}{y}_3^n+c_{\alpha z}{y}_4^n{)}^2\\ & +\frac{w_{\alpha }}{2(y_1^n{)}^2{c}_s^2}((y_2^n{)}^2+(y_3^n{)}^2+(y_4^n{)}^2)]\end{array}\text{\hspace{1em}(17)}$$

格子张量的一些性质²
$$\begin{array}{rlr}& \sum _{\alpha }{w}_{\alpha }=1.& \alpha \in [0,1,...,q-1]\\ & \sum _{\alpha }{w}_{\alpha }{c}_{\alpha \beta }=0.& \beta \in [x,y,z]\\ & \sum _{\alpha }{w}_{\alpha }{c}_{\alpha \beta }{c}_{\alpha \gamma }=c_s^2{\delta }_{\beta \gamma }.& \beta ,\gamma \in [x,y,z]\end{array}$$
注: $q$ 为离散速度数目

基于格子张量的性质，式 (17) 右侧仅剩下 $\sum _{\alpha }{w}_{\alpha }$ ，即：

$${\left[\frac{\partial y_1^{\ast }}{\partial y_1}\right]}^n=1$$

同理，其他项也可进行类似分析，并得出

$${\left[\frac{\partial y_{\lambda }^{\ast }}{\partial y_{\kappa }}\right]}^n={\delta }_{\lambda \kappa }$$

其中 ${\delta }_{\lambda \kappa }$ 为狄拉克函数，且 $\lambda ,\kappa \in [1,2,3,4]$ 。所以我们可以得到 ${\overline{\mathrm{d}\mathbf{Y}}}^{\ast }={\overline{\mathrm{d}\mathbf{Y}}}^n$ 。

校正步

令 ${\mathrm{d}{y}_i}^{n+1}={\overline{\mathrm{d}{y}_i}}^{n+1}\exp (j\mathbf{k}\cdot \mathbf{r})$（$i\in [1,2,3,4]$）。由于校正步（式 (9-2) ）中并未对 $\rho$ 做修改，易得：

$${\overline{\mathrm{d}{y}_1}}^{n+1}={\overline{\mathrm{d}{y}_1}}^{\ast }={\overline{\mathrm{d}{y}_1}}^n$$

根据式 (9-2) 的校正步，以及 $f_{\alpha }^{(neq)}(\mathbf{r}+\frac{\Delta t}{2}{\mathbf{c}}_{\alpha },t-\frac{\Delta t}{2})$ 的差分表示，有

$$\begin{array}{r}{y}_2^{n+1}=y_2^{\ast }-(\tau -1)y_2^n+(\tau -1)\sum _{\alpha }{c}_{\alpha x}{f}_{\alpha }^{(eq)}(\mathbf{r}+{\mathbf{c}}_{\alpha }\Delta t,t)\\ y_3^{n+1}=y_3^{\ast }-(\tau -1)y_3^n+(\tau -1)\sum _{\alpha }{c}_{\alpha y}{f}_{\alpha }^{(eq)}(\mathbf{r}+{\mathbf{c}}_{\alpha }\Delta t,t)\\ y_4^{n+1}=y_4^{\ast }-(\tau -1)y_4^n+(\tau -1)\sum _{\alpha }{c}_{\alpha z}{f}_{\alpha }^{(eq)}(\mathbf{r}+{\mathbf{c}}_{\alpha }\Delta t,t)\end{array}\text{\hspace{1em}(18)}$$

这里同样以 $y_2^{n+1}$ 的计算为例，记

$$\begin{array}{rl}S& =\sum _{\alpha }{c}_{\alpha x}{f}_{\alpha }^{(eq)}(\mathbf{r}+{\mathbf{c}}_{\alpha }\Delta t,t)\\ & =\sum _{\alpha }{c}_{\alpha x}{w}_{\alpha }{y}_{1(-\alpha )}^{\ast }\\ & +\sum _{\alpha }\frac{c_{\alpha x}{w}_{\alpha }}{c_s^2}(c_{\alpha x}{y}_{2(-\alpha )}^{\ast }+c_{\alpha y}{y}_{3(-\alpha )}^{\ast }+c_{\alpha z}{y}_{4(-\alpha )}^{\ast })\\ & +\sum _{\alpha }\frac{c_{\alpha x}{w}_{\alpha }}{2c_s^4{y}_{1(-\alpha )}^{\ast }}(c_{\alpha x}{y}_{2(-\alpha )}^{\ast }+c_{\alpha y}{y}_{3(-\alpha )}^{\ast }+c_{\alpha z}{y}_{4(-\alpha )}^{\ast }{)}^2\\ & -\sum _{\alpha }\frac{c_{\alpha x}{w}_{\alpha }}{2c_s^2{y}_{1(-\alpha )}^{\ast }}((y_{2(-\alpha )}^{\ast }{)}^2+(y_{3(-\alpha )}^{\ast }{)}^2+(y_{4(-\alpha )}^{\ast }{)}^2{)}^2\end{array}$$

其中下标 $(-\alpha )$ 代表来自 $\mathbf{r}+{\mathbf{c}}_{\alpha }\Delta t$ 的值。则

$$y_2^{n+1}=y_2^{\ast }+(\tau -1)(S-y_2^n)\text{\hspace{1em}(19)}$$

所以其差分表示为

$$\mathrm{d}{y}_2^{n+1}=\mathrm{d}{y}_2^{\ast }+(\tau -1)(\mathrm{d}S-\mathrm{d}{y}_2^n)\text{\hspace{1em}(20)}$$

对于 $\mathrm{d}S$，有

$$\mathrm{d}S={\left(\frac{\partial S}{\partial y_1}\right)}^{\ast }\mathrm{d}{y}_1^{\ast }+{\left(\frac{\partial S}{\partial y_2}\right)}^{\ast }\mathrm{d}{y}_2^{\ast }+{\left(\frac{\partial S}{\partial y_3}\right)}^{\ast }\mathrm{d}{y}_3^{\ast }+{\left(\frac{\partial S}{\partial y_4}\right)}^{\ast }\mathrm{d}{y}_4^{\ast }$$

以 ${\left(\frac{\partial S}{\partial y_1}\right)}^{\ast }$ 为例，有

$$\begin{array}{rl}{\left(\frac{\partial S}{\partial y_1}\right)}^{\ast }& =\sum _{\alpha }{c}_{\alpha x}{w}_{\alpha }\\ & -\sum _{\alpha }\frac{c_{\alpha x}{w}_{\alpha }}{2c_s^4(y_{1(-\alpha )}^{\ast }{)}^2}(c_{\alpha x}{y}_{2(-\alpha )}^{\ast }+c_{\alpha y}{y}_{3(-\alpha )}^{\ast }+c_{\alpha z}{y}_{4(-\alpha )}^{\ast }{)}^2\\ & +\sum _{\alpha }\frac{c_{\alpha x}{w}_{\alpha }}{2c_s^2(y_{1(-\alpha )}^{\ast }{)}^2}((y_{2(-\alpha )}^{\ast }{)}^2+(y_{3(-\alpha )}^{\ast }{)}^2+(y_{4(-\alpha )}^{\ast }{)}^2{)}^2\end{array}$$

与预报步的做法类似，将 $y_{i(-\alpha )}^{\ast }$ 近似为 $y_i^{\ast }$，得到：

$$\begin{array}{rl}{\left(\frac{\partial S}{\partial y_1}\right)}^{\ast }& =\sum _{\alpha }{c}_{\alpha x}{w}_{\alpha }\\ & -\sum _{\alpha }\frac{c_{\alpha x}{w}_{\alpha }}{2c_s^4(y_1^{\ast }{)}^2}(c_{\alpha x}{y}_2^{\ast }+c_{\alpha y}{y}_3^{\ast }+c_{\alpha z}{y}_4^{\ast }{)}^2\\ & +\sum _{\alpha }\frac{c_{\alpha x}{w}_{\alpha }}{2c_s^2(y_1^{\ast }{)}^2}((y_2^{\ast }{)}^2+(y_3^{\ast }{)}^2+(y_4^{\ast }{)}^2{)}^2\end{array}$$

此时所有 $y_i^{\ast }$ 均可被视为常数提取出来。根据格子张量的性质，可得： ${\left(\frac{\partial S}{\partial y_1}\right)}^{\ast }=0$。

同理
$${\left(\frac{\partial S}{\partial y_2}\right)}^{\ast }=1,{\left(\frac{\partial S}{\partial y_3}\right)}^{\ast }={\left(\frac{\partial S}{\partial y_4}\right)}^{\ast }=0$$

此时 $\overline{\mathrm{d}S}={\overline{\mathrm{d}{y}_2}}^{\ast }$ ，则式 (20) 可化简为：

$${\overline{\mathrm{d}{y}_2}}^{n+1}={\overline{\mathrm{d}{y}_2}}^{\ast }={\overline{\mathrm{d}{y}_2}}^n$$

同理，对 $y_3^{n+1}$ 和 $y_4^{n+1}$ 计算可得： ${\overline{\mathrm{d}{y}_3}}^{n+1}={\overline{\mathrm{d}{y}_3}}^n$ 和 ${\overline{\mathrm{d}{y}_4}}^{n+1}={\overline{\mathrm{d}{y}_4}}^n$ 。

记 $\mathrm{d}{\mathbf{Y}}^{\mathrm{n}}=(\mathrm{d}{y}_1^n,\mathrm{d}{y}_2^n,\mathrm{d}{y}_3^n,\mathrm{d}{y}_4^n{)}^T$， $\mathrm{d}{\mathbf{Y}}^{\mathrm{n}+1}=(\mathrm{d}{y}_1^{n+1},\mathrm{d}{y}_2^{n+1},\mathrm{d}{y}_3^{n+1},\mathrm{d}{y}_4^{n+1}{)}^T$ ，那么：

$${\overline{\mathrm{d}\mathbf{Y}}}^{n+1}=\mathbf{I}{\overline{\mathrm{d}\mathbf{Y}}}^n$$

其中 $\mathbf{I}$ 为单位矩阵。显然，特征矩阵的特征值全是 1 ，代表着该格式无条件稳定。

---

1. Chen, Z., Shu, C., Wang, Y., Yang, L. M., & Tan, D. (2017). A Simplified Lattice Boltzmann Method without Evolution of Distribution Function. Advances in Applied Mathematics and Mechanics, 9(1), 1–22. DOI:10.4208/aamm.OA-2016-0029 ↩︎

2. Z. Chen, C. Shu, D. Tan; Three-dimensional simplified and unconditionally stable lattice Boltzmann method for incompressible isothermal and thermal flows. Physics of Fluids. 1 May 2017; 29 (5): 053601. DOI:10.1063/1.4983339 ↩︎ ↩︎ ↩︎

3. Chen Z, Shu C, Tan D, Wu C. On improvements of simplified and highly stable lattice Boltzmann method: Formulations, boundary treatment, and stability analysis. Int J Numer Meth Fluids. 2018; 87: 161–179. DOI:10.1002/fld.4485 ↩︎

4. Chen, Z., & Shu, C. (2020). Simplified and Highly Stable Lattice Boltzmann Method. WORLD SCIENTIFIC. DOI:10.1142/12047 ↩︎

5. 本文中所有 $\nabla$ 均表示空间偏导 $\frac{\partial }{\partial \mathbf{r}}$。 ↩︎

6. C. Shu, Y. Wang, C. Teo, and J. Wu. Development of lattice Boltzmann flux solver for simulation of incompressible flows. Advances In Applied Mathematics And Mechanics. 6, 436 (2014). ↩︎