格子Boltzmann方法（LBM）基础理论推导以及与Taichi联合的原因

[Tau]

LBM理论基础推导以及与Taichi联合的原因

• 笔者身份：LBM初学者，并行计算纯小白，某不知名高校研一学生在读

• 定位：我认为与其认为这是一次分享，不如认为它是我的一个学习笔记

• 目的：我希望可以通过这次分享，使得许多初学LBM的同学可以简单了解其理论基础，并且了解我为何在对比之后选择使用Taichi语言。当然由于本人也只是学生，因此能力有限，文中无可避免会存在一些纰漏，恳请大家批评指正。

• 联系方式：有问题可以随时评论或者直接邮件交流，邮箱：Lizt1191@gmail.com

• 参考资料：

  [1] 郭照立, 郑楚光. 格子 Boltzmann 方法的原理及应用[M]. 科学出版社, 2009.

  [2] 何雅玲, 王勇, 李庆. 格子 Boltzmann 方法的理论及应用[M]. 科学出版社, 2009.

  [3] Guo Z, Shu C. Lattice Boltzmann method and its application in engineering[M]. World Scientific, 2013.

  ---

一， 流体力学基本方程

1. 雷诺输运方程

 众所周知，描述流体运动过程一般有两种法：拉格朗日的法和欧拉法，前者的主要思路是跟随流体质点去研究流体运动，而后者则是着眼于空间坐标去研究流体运动，两种方法的差异于联系在任意流体力学相关的书籍中应该都有所描述，此处就不再赘述。而在推导流体力学基本方程的过程中，为了将针对于系统的物理量转换到易于研究的欧拉参考系当中，我们需要引入雷诺输运方程，其作用是将一个流体系统的拉格朗日变化率表示为欧拉导数。其具体推导过程可参见前文所列的参考文献，此处只给出最终结果作为前备知识。

\frac{\mathrm{D}}{\mathrm{D} t} \int_{V(t)} \phi \mathrm{d} V=\int_{V(t)} \frac{\partial \phi}{\partial t} \mathrm{~d} V+\int_{A(t)} \phi u \cdot n \mathrm{~d}

2. 连续性方程

 由质量守恒： \frac{\mathrm{D} m}{\mathrm{D} t}=\frac{\mathrm{D}}{\mathrm{D} t} \int_{V(t)} \rho \mathrm{d} V=0 ，又由雷诺输运方程

\frac{\mathrm{D}}{\mathrm{D} t} \int_{V(t)} \rho d V=\int_{V(i)}\left[\frac{\partial \rho}{\partial t}+\nabla \cdot(\rho u)\right] \mathrm{d} V=0

\Downarrow

\frac{\partial \rho}{\partial t}+\nabla \cdot(\rho u)=0 \Leftrightarrow \frac{D \rho}{D t}+\rho \nabla \cdot u=0

3. 动量方程

 由牛顿第二定律，系统内动量变化率=作用在系统的总质量力+\Sigma 作用在系统上的表面力

\Downarrow

\frac{D}{D t} \int_{V(t)}(\rho u) d V=\int_{V(t)} \rho \boldsymbol{F} d V+\int_{A(t)} \sigma \cdot n d A

\Downarrow

\int_{V(t)} \rho \frac{D u t}{D t} d V=\int_{V(t)} \rho \boldsymbol{F} d V+\int_{V(t)} \nabla \cdot \sigma d V

\Updownarrow

\int_{V(t)}\left(\rho \frac{D u t}{D t}-\rho \boldsymbol{F}-\nabla \cdot \sigma\right) d V=0

\Updownarrow

\rho \frac{D u}{D t}=\rho \boldsymbol{F}+\nabla \cdot \sigma

又 \frac{D \rho}{D t}+\rho \nabla \cdot u=0

\Downarrow

\frac{\partial(\rho u)}{\partial t}+\nabla \cdot(\rho d d)=\rho \boldsymbol{F}+\nabla \cdot \sigma

其中 \boldsymbol{F} 为质量力， \sigma 为应力张量，对于牛顿流体，由斯托克斯的三项假设，可将 \sigma 写作 \sigma=-p\boldsymbol{I}+\boldsymbol{\Pi} ，其中 \boldsymbol{I} 为单位二阶张量， \boldsymbol{\Pi}=\mu\left[\left(\nabla u+(\nabla u)^{\mathrm{T}}\right)-\frac{2}{3}(\nabla \cdot u) \boldsymbol{I}\right] ，该公式可适用于除高温和高频声波等极端情况外的一般气体运动过程中，具体推导过程详见参考书籍。将应力张量的表达式代入可得

\rho \frac{\mathrm{D} u}{\mathrm{D} t}=\rho \boldsymbol{F}-\nabla p+\nabla \cdot\left[\mu\left(\nabla u+(\nabla u)^{\mathrm{T}}\right)\right]+\nabla\left(\mu^{\prime} \nabla \cdot u\right)

其中 \mu' 为第二粘度系数， \mu'=-\frac{2}{3} \mu +\mu_B ， \mu_B 为第二粘性系数，一般取0，当温度变化很小的时候，可以认为 \mu 和 \mu' 均为常数 ，故有：

对于可压缩流体：\frac{\partial(\rho u)}{\partial t}+\nabla \cdot(\rho u u)=\rho \boldsymbol{F}-\nabla p+\mu \nabla^2 u+\left(\mu+\mu^{\prime}\right) \nabla(\nabla \cdot u)

对于不可压流体：\frac{\partial u}{\partial t}+\nabla \cdot(u u)=\boldsymbol{F}-\frac{1}{\rho} \nabla p+\nu \nabla^2 u

二，Boltzmann方程

 Boltzmann 方程是统计力学中用以描述非平衡态分布函数演化规律的方程，从统计力学的角度来讲，在任何一个宏观系统中，其中每个分子的微观运动都遵守力学规律，因此只需要算出大量粒子的个别运动,就可以确定系统的宏观参数；从另–个方面来考虑,我们可以不去确定每个分子的运动状态，而是求出每一个分子处在某一状态下的概率，之后通过统计方法得出系统的宏观参数，这是Boltzmann方程的基本思想。

 我们首先需要定义分布函数 f ，其是空间位置矢量 r(x, y, z) 分子速度矢量 \xi\left(\xi_x, \xi_y, \xi_z\right) 及时间 t 的函数。任一分子假如在一定时间内无碰撞，则有

f(r+\mathrm{d} r, \xi+a \mathrm{~d} t, t+\mathrm{d} t) \mathrm{d} r \mathrm{~d} \xi=f(r, \xi, t) \mathrm{d} r \mathrm{~d} \xi

由泰勒展开 \Rightarrow \frac{\partial f}{\partial t}+\xi \cdot \frac{\partial f}{\partial r}+a \cdot \frac{\partial f}{\partial \xi}=0

\Updownarrow

\left(\frac{\partial f}{\partial t}\right)_{\text {运动 }}=-\xi \cdot \frac{\partial f}{\partial r}-\boldsymbol{a} \cdot \frac{\partial f}{\partial \xi}

又\left(\frac{\partial f}{\partial t}\right)=\left(\frac{\partial f}{\partial t}\right)_{\text {运动 }}+\left(\frac{\partial f}{\partial t}\right)_{\text {碰撞 }}

 由钢球碰撞模型， \left(\frac{\partial f}{\partial t}\right)_{\text {碰撞 }}=\iint\left(f^{\prime} f_1^{\prime}-f f_1\right) d_{\mathrm{D}}^2|g| \cos \theta \mathrm{d} \Omega \mathrm{d} \xi_1 ，其中 \int \mathrm{d} \Omega=\int_0^{2 \pi} \mathrm{d} \phi \int_0^{\pi / 2} \sin \theta \mathrm{d} \theta ， \int \mathrm{d} \xi_1=\iiint_{-\infty}^{+\infty} \mathrm{d} \xi_{1 x} \mathrm{~d} \xi_{1 y} \mathrm{~d} \xi_{1 z} 。将其带入可得：

\frac{\partial f}{\partial t}+\xi \cdot \frac{\partial f}{\partial r}+a \cdot \frac{\partial f}{\partial \xi}=\iint\left(f^{\prime} f_1^{\prime}-f f_1\right) d_{\mathrm{D}}^2|g| \cos \theta \mathrm{d} \Omega \mathrm{d} \xi_1

即为一般形式的Boltzmann方程，是气体动理学基本方程，等号右端的项称为碰撞项，一般用 J\left(f f_1\right) 表示。一般情况下要获得Boltzmann方程的解是不可能的，不过单组分单原子气体在不受外力作用的情况下，可以得到他的一个理论解：Maxwell分布，即为单组分单原子气体的平衡态分布，具体形式如下：

f=n \frac{1}{\left(2 \pi R_{\mathrm{g}} T\right)^{D / 2}} e ^{-\frac{(\xi-u)^2}{2 R_{\mathrm{g}} T}}

 其整体的推导思路：首先从平衡态时分子碰撞不变出发，由碰撞不变量的线性组合是仅有的求和不变量，可以设出平衡分布函数，后用分子数密度和温度的表达式可以推导出平衡分布函数的系数，进而得到Maxwell分布，具体过程与相关参数详见参考书。

三，两者之间的联系与推导

 首先提出Boltzmann方程的矩方程的概念：将一个与分子相关联的量 \phi （ \phi 是 r、g和t 的一个任意函数，其可以是气体分子的质量、动量、能量等）与 Boltzmann方程相乘并逐项对速度空间积分，即可得到 Boltzmann方程的矩,这称之为Boltzmann方程的矩方程，也将其称之为 Maxwell输运方程：

\int \phi\left(\frac{\partial f}{\partial t}+\xi \cdot \frac{\partial f}{\partial r}+a \cdot \frac{\partial f}{\partial \xi}\right) \mathrm{d} \xi=\int \phi J\left(f f_1\right) \mathrm{d} \xi

\Downarrow

\frac{\partial(n \bar{\phi})}{\partial t}+\frac{\partial}{\partial r} \cdot(n \overline{\phi\xi})-n\left(\frac{\partial \bar{\phi}}{\partial t}+\overline{\xi \cdot \frac{\partial \phi}{\partial r}}+a \cdot \overline{\partial \phi}\right)=\int \phi J\left(f f_1\right) d \xi

其中若 \phi 只和分子速度$\xi$ 有关而与位置矢量 r 与时间 t 都无关，则 \frac{\partial \bar{\phi}}{\partial t}+\overline{\xi \cdot \frac{\partial \phi}{\partial r}} 就可以忽略，若某物理量在碰撞前后守恒，则 \int \phi J\left(f f_1\right) \mathrm{d} \xi=0 ，而很显然，质量，动量，能量，这几个量为碰撞不变量，因此：

• 若 \phi=m ，代入得 :

  质量\longrightarrow质量守恒方程

\frac{\partial \rho}{\partial t}+\nabla \cdot(\rho u)=0

其中 \rho=nm 为气体密度， u=\frac{1}{n} \int \xi f d \xi 为气体宏观平均速度

• 若 \phi=m\xi ，代入得：

  动量\longrightarrow动量守恒方程

\frac{\partial}{\partial t}(\rho u)+\frac{\partial}{\partial r} \cdot \int m \xi \xi f \mathrm{~d} \xi-\rho a=0

其中\int m \xi \xi f d \xi=-\sigma+\rho uu

\Downarrow

\frac{\partial}{\partial t}(\rho u)+\nabla \cdot(\rho u t)=\nabla \cdot \sigma+\rho a

• 若 \phi=\frac{1}{2}m\xi^2 ，代入得：

  能量\longrightarrow能量守恒方程

\frac{\partial}{\partial t} \int \frac{1}{2} m \xi^2 f \mathrm{~d} \xi+\frac{\partial}{\partial r} \cdot \int \frac{1}{2} m \xi^2 \xi f \mathrm{~d} \xi-a \cdot \int f m \xi \mathrm{d} \xi=0

其中\int \frac{1}{2} m \xi^2 f \mathrm{~d} \xi =\rho \varepsilon+\frac{1}{2} \rho u^2 ，且内能\varepsilon=\frac{1}{2 n} \int c^2 f \mathrm{~d} \xi

\int \frac{1}{2} m \xi^2 \xi f d \xi=q+\rho\left(\varepsilon+\frac{1}{2} u^2\right) u+\sigma \cdot u，且热流矢量q=\int \frac{1}{2} m c^2 c f d \xi

\Downarrow

\frac{\partial}{\partial i}\left(\rho \varepsilon+\frac{1}{2} \rho u^2\right)+\nabla \cdot\left[\left(\rho \varepsilon+\frac{1}{2} \rho u^2\right) u\right]=-\nabla \cdot q+\nabla \cdot(\sigma \cdot u)+\rho a

四，两者的差异以及各自的应用范围

 当研究对象的尺度远远大于粒子的结构尺度的时候，就可以认为连续介质假设成立，因此其流体力学领域被广泛的采用并取得了极大的成功。例如飞机在空气中的运动，轮船在水中的运动，由于其特征长度远远大于粒子的结构尺度，就可以空气和水在这个系统中就可以被认为是连续介质；然而,在研究高空稀薄气体中的物体运动时,因分子平均自由程很大，与物体特征长度尺度相比为同阶量，这时便不能视高空稀薄气体为连续介质；同样，血液在微血管中的运动也不能被当作连续介质等等。这些都是连续介质应用范围的限制。

 对于连续介质假设不能适用或宏观方程难以描述的系统，采用微观的分子动力学或介观的气体动理论(gas kinetic theory)进行描述更为恰当和可行。但是，这两类方法又会带来其他问题。分子动力学方法需要跟踪多个粒子的运动，计算量非常之大，对计算机的存储量和计算速度都有很高的要求。而气体动理论的基本方程正是前文所推导的连续的Boltzmann方程，这是个比N-S方程组更加复杂的微分积分方程，其准确求解更加不现实，因此需要对其加以离散化来构造微观或者介观的控制方程。

 这些离散格式都是基于这样的事实,即流体的宏观运动是流体分子微观热运动的统计平均结果,宏观行为对每个具体分子的运动细节是不敏感的。Navier-Stokes方程组所描述的守恒定律与微观粒子的运动规律是一致的,流体分子内部的相互作用的差别反映在Navier-Stokes方程组的输运系数上。因此,人们可以构造微观或介观的模型,使之在遵循基本守恒律的前提下尽可能地简洁，也因此进一步的出现了各种的离散格式

 Boltzmann方程在整个稀薄气体动力学中占据中心地位：在自由分子流区中,通常应用无碰撞项的Boltzmann方程；在滑移区，则应用从解 Boltzmann方程的Chapman-Enskog展开所得到的一阶近似Navier-Stokes方程组或二阶近似 Burnett方程组；在过渡区则利用 Boltzmann方程或用与它等价的方法求解气体流动的问题。

五，由Boltzmann方程到格子Boltzmann方程（LBE）

 由前文所述，我们对各个公式的关系做一总结：Boltzmann方程是气体动理论的基本方程，也是整个稀薄气体动力学的基本方程；Maxwell分布是Boltzmann方程在单组分单原子气体时的平衡态分布；而由于Boltzmann方程是一个复杂的微分积分方程，它的直接求解通常是不可能的，于是人们便提出了许多近似解法,例如用 Boltzmann-BGK方程代替 Boltzmann方程，其主要针对于Boltzmann方程的碰撞项利用一个简单的碰撞算子 \Omega_f 进行了近似处理，使得其可以求解，本节将主要围绕此部分展开。

 最简单的一个碰撞算子 \Omega_f 认为碰撞的效应是改变 f 使其趋于Maxwell平衡态分布 f^{eq} 。设改变率的大小和 f^{eq} 与 f 的差值成正比,比例系数为 \nu 。 \nu 是一个与分子速度 \xi 无关的常数，于是碰撞算子 \Omega_f 便可写作：

\Omega_f=v\left[f^{\mathrm{eq}}(r, \xi)-f(r, \xi, t)\right]

\Downarrow

Boltzmann-BGK方程：\frac{\partial f}{\partial t}+\xi \cdot \frac{\partial f}{\partial r}+a \cdot \frac{\partial f}{\partial \xi}=v\left(f^{\text {eq }}-f\right)

 BGK近似，使得 Boltzmann方程得以线性化，大大简化了方程的求解。如果引人碰撞时间 \tau_0=\frac{1}{\nu} ，其物理意义为每两次碰撞之间的时间间隔，也称松弛时间或者弛豫时间，在物理学中，将非平衡态分布向平衡态分布趋近的过程称为弛豫过程。因此BGK近似是把气体从当前状态向平衡态的过渡看作是个简单的弛豫过程，因其简单而被广泛应用。但它终究是用近似的方式替代了以坚实物理现实为基础的准确项，因而必然存在有缺陷，不过对于偏离平衡态不远的小扰动问题，BGK 近似的精度可以满足要求。在其分布函数的两边同时乘以 m 来使其表示粒子的密度分布，用松弛时间可以将其改写为：

\frac{\partial(m f)}{\partial t}+\xi \cdot \nabla(m f)+a \cdot \nabla_{\xi}(m f)=-\frac{1}{\tau}\left(m f-f^{e q}\right)

其对应的Maxwell平衡态分布 f^{eq} 变为：

f^{e q}=\rho \frac{1}{\left(2 \pi R_g T\right)^{D / 2}} \exp \left[-\frac{(\xi-u)^2}{2 R_g T}\right]

 对Boltzmann方程进行BGK近似之后使得其线性化，使得方程可以求解，接下来为了数值求解Boltzmann-BGK方程，需要对其进行离散化，将其速度离散后可得：

\frac{\partial f_\alpha}{\partial t}+e_\alpha \cdot \nabla f_\alpha=-\frac{1}{\tau_0}\left(f_\alpha-f_\alpha^{\mathrm{eq}}\right)+F_\alpha

其中 F_\alpha 为外力项， F_\alpha=\left(a \cdot \nabla_\xi f\right)_\alpha

 当流速较低时，局部的平衡态分布函数 f_\alpha^{eq} 可以表示为：

f_\alpha^{\mathrm{eq}}=\rho \omega_a\left[1+\frac{e_\alpha \cdot u}{R_{\mathrm{g}} T}+\frac{\left(e_a \cdot u\right)^2}{2 R_{\mathrm{g}}^2 T^2}-\frac{u^2}{2 R_{\mathrm{g}} T}\right]+O\left(u^3\right)

其中 \omega_\alpha=\left(2 \pi R_{\mathrm{g}} T\right)^{-D / 2} \exp \left[-e_\alpha^2 /\left(2 R_{\mathrm{g}} T\right)\right] ，在对其在时间和空间上离散之后得到完全离散化的格子Boltzmann-BGK方程：

f_\alpha\left(r+e_\alpha \delta_t, t+\delta_t\right)-f_\alpha(r, t)=-\frac{1}{\tau}\left[f_\alpha(r, t)-f_\alpha^{\mathrm{eq}}(r, t)\right]+\delta_t F_\alpha(r, t)

由参考书中的推导可知，其基本可以达到二阶数值计算精度，并且时间和空间的离散并不独立，而是由粒子的离散速度来将其联系起来，即 \delta_x=e_{\alpha x} \delta_t ，这一特性为我们在物理空间将粒子运动分为迁移和碰撞两个相对过程提供了条件，即粒子在两个时间步之间由一个节点恰好运动到对应的相邻节点，并在该节点上与其他粒子碰撞,这使格子Boltzmann方法具备了很好的并行特性和较强复杂边界的处理能力，也正是这一点，为LBM和Taichi的联用创造了极大的优势。

 对于一个完整的格子Boltzmann模型来说，其一般会由三部分构成：离散速度模型，平衡态分布函数，分布函数的演化方程。构造格子Boltzmann模型的关键是选择合适的平衡态分布函数，而平衡态分布函数的具体形式又与离散速度模型的构造有关，离散速度的对称性决定了相应的格子Boltzmann模型能否还原到所要求解的宏观方程。因此,离散速度模型的构造显得至关重要：如果离散速度的个数太少可能导致某些应当守恒的物理量不满足守恒律，太多则可能造成计算上的浪费。由Qian等人提出的DdQm（d维空间，m个离散速度）模型是格子Boltzmann方法的基本模型，因此在此简单列举最常用的D2Q9模型的完整描述。

 演化方程：

f_\alpha\left(r+e_\alpha \delta_t, t+\delta_t\right)-f_\alpha(r, t)=-\frac{1}{\tau}\left[f_\alpha(r, t)-f_\alpha^{\mathrm{eq}}(r, t)\right]+\delta_t F_\alpha(r, t)

 其离散速度配置为：

\mathbb{e}_\alpha= \begin{cases}(0,0) & \alpha=0 \\ c\left(\cos \left[(\alpha-1) \frac{\pi}{2}\right], \sin \left[(\alpha-1) \frac{\pi}{2}\right]\right) & \alpha=1,2,3,4 \\ \sqrt{2} c\left(\cos \left[(2 \alpha-1) \frac{\pi}{4}\right], \sin \left[(2 \alpha-1) \frac{\pi}{4}\right]\right) & \alpha=5,6,7,8\end{cases}

其中 c=\delta_x / \delta_t

 平衡态分布函数可写作如下形式：

f_\alpha^{\mathrm{eq}}=\rho w_a\left[1+\frac{e_\alpha \cdot u}{c_{\mathrm{s}}^2}+\frac{\left(e_\alpha \cdot u\right)^2}{2 c_{\mathrm{s}}^4}-\frac{u^2}{2 c_{\mathrm{s}}^2}\right]

其中 C_s 为格子声速， \omega_\alpha 为权系数

 若要恢复宏观方程，则平衡态分布函数需满足以下矩方程:

\begin{gathered} \sum_\alpha f_\alpha^{\mathrm{eq}}=\rho \\ \sum_\alpha f_a^{\mathrm{eq}} \mathbb{e}_\alpha=\rho u\\ \sum_\alpha f_\alpha^{\mathrm{eq}} e_{\alpha i} e_{\alpha j}=\rho u_i u_j+p \delta_{i j} \end{gathered}

分别代入整理可得权系数和格子声速的计算式：

w_0=\frac{4}{9}, \quad w_1=\frac{1}{9}, \quad w_2=\frac{1}{36}, \quad c_{\mathrm{s}}^2=\frac{c^2}{3}

并有： p=\rho c_s^2 ，模型的宏观密度和速度如下定义：

\begin{gathered} \rho=\sum_\alpha f_\alpha \\ u=\frac{1}{\rho} \sum_\alpha f_\alpha e_\alpha \end{gathered}

 通过C-E展开，若令运动粘度系数 \nu=c_{\mathrm{s}}^2\left(t-\frac{1}{2}\right) \delta_t ，则可将上述的格子Boltzmann方程在不考虑外力的情况下得到其对应的宏观方程：

\frac{\partial \rho}{\partial t}+\nabla \cdot(\rho \boldsymbol{u})=0

\frac{\partial(\rho u)}{\partial t}+\nabla \cdot(\rho u u)=-\nabla p+\nabla \cdot\left[\rho \nu\left(\nabla u+(\nabla u)^{\mathrm{T}}\right)-\frac{\nu}{c_{\mathrm{s}}^2} \nabla \cdot(\rho u u u)\right]

 与标准的可压缩N-S方程组相比，虽然连续性方程完全一致，但动量方程存在偏差，此时若流体的密度恒定为常数，且对于低马赫数流动来说， \nabla \cdot \text { (puuu) } 这一项也就可以忽略，此时其与N-S方程组可以保持一致

 需要注意的是，上述推导出来的格子Boltzmann基本模型所导出的宏观方程从形式上来说是可压缩的N-S方程组，当用其模拟不可压缩流动时，会存在一定的压缩性误差，为降低和消除这些误差，需要对其进行一定的改进，使得其恢复的宏观方程与不可压的N-S方程组相同或非常接近，后面也将会根据本人的应用来对其进行展示。

六，LBM常见的边界格式

 关于LBM的边界格式这块，限于本人的能力，此处只简要介绍几种我将会使用到的边界格式，且不进行理论的推导，只进行公式的展示：

 启发式格式是一种根据流场的流动特点对边界上的格点直接设定分布函数的边界条件，一些有比较明显宏观特性的边界就常常采用启发式格式。例如较为常见的周期性流场，其假设当一个流体点从流场一边的边界离开流场的时候，在下一时间步就会从流场的另外一边再流入流场，很显然这种情况下流场的质量和动量是守恒的，这种边界格式常常使用于一些可以认为在某个方向无穷大的流动。假设其在x方向的流动的周期是，则可以将该格式表示为：

f_i\left(0, y, z, t+\delta_t\right)=f_i^{\prime}\left(L_x, y, z, t\right), \quad f_i\left(L_x, y, z, t+\delta_t\right)=f_i^{\prime}(0, y, z, t)

 对于无滑移壁面的边界的处理常常采用启发式格式中的Half-Way反弹格式（半反弹格式），其是在标准反弹格式的基础之上进行改进得到的，虽然和标准反弹格式一样，都是在与固体边界相邻的流体格点上执行碰撞，但是由于半反弹格式的固体边界是放置在流体格点与边界固体格点的中间处，即 \left(x_f+x_w\right) / 2 处，当在流体格点（ x_f ）处进行碰撞之后，按照标准的碰撞-迁移过程向壁面迁移的粒子经过 \delta_t/2 之后到达壁面，与壁面进行碰撞并反转其迁移路线，再经 \delta_t/2 后到达格点 x_f 处，因此必有：

f_{\bar{\imath}}^{\prime}\left(x_b, t\right)=f_i^{\prime}\left(x_f, t\right)

其中 e_{\bar{l}}=-e_i （ e_i 指向壁面）， x_b 是壁面格点， x_f=x_b-c_i\delta_t 是流体格点。同时半反弹边界格式具有二阶精度，而标准反弹只有一阶精度，因此在之后的代码中主要应用半反弹的边界格式。

 由上述碰撞过程的描述，我们可以发现半反弹格式具有清晰的物理图像，且不增加任何额外的计算量，对不规则的复杂边界的处理尤为适用。反弹格式作为LBE方法的一个显著的特点，使得其在模拟复杂流-固相互作用的流动（如多孔介质，变形边界等）具有很大的优势。

 在启发式格式和动力学格式之外还有一种外推格式，与前两种格式不同的是，外推格式主要是从求解数学物理方程的角度来发展并推导出来的，因此其较为适用于一些非常规的边界条件（如包含梯度信息的边界）。受Chen的外推方法和Zou的非平衡态反弹方法的启发，吸收两种方法的优点，Guo等于2002年提出一种新的边界处理方法，即非平衡态外推格式（Nonequilibrium Extrapolation Scheme）。其基本思想是：将边界节点上的分布函数分解为平衡态和非平衡态两部分，平衡态部分由具体的宏观边界条件来构造新的平衡态分布函数来近似获得，非平衡态部分则使用一阶精度的外推确定，具体形式可以写作：

 我们用速度边界来说明该格式的基本原理，2，5，6方向指向流体侧，0，1，3方向处于边界上，4，7，8方向则指向流场外部，也就是固体侧。在边界格点上的2、5、6方向的分布函数无法通过碰撞与迁移的方式获得，因此需要外推出其此时的分布函数，其计算格式如下：

f_i\left(\boldsymbol{x}_b, t\right)=f_i^{e q}\left(\rho\left(\boldsymbol{x}_f, t\right), \boldsymbol{u}_w\right)+\left[f_i\left(\boldsymbol{x}_f, t\right)-f_i^{e q}\left(\boldsymbol{x}_f, t\right)\right], \quad i=2,5,6

 郭照立等证明，非平衡态外推格式整体具有二阶精度，和LBGK模型的精度一致。并且在具有较好数值稳定性的基础之上，还有适用范围广，容易实现，计算简单等优点，因此得到了广泛的采用。

七，与Taichi联用的原因

 在前面几节中已经大致介绍了LBM的一些独特的优点以及其适用的范围，并且由于其采用的碰撞迁移规则来求解离散后的LB-BGK方程，从而避免了复杂的矩阵方程组的求解，为采用并行计算带来了非常大的便利，很多学者就采用OpenMP，MPI和CUDA编程以实现高性能并行计算，但是由于其本身的复杂性，再加上学习其所需要的基础知识多且庞杂，一不小心就会因为一步的疏忽而满盘报错，而对于非高性能计算相关从业者而言，这些知识并非是我们所关心的内容，我只是需要一个简单易用的工具来帮我实现并行计算这一步，相比较最大化的效率而言，其简洁性和所要花费的沉没成本则是我更加关心的点。我本来这段时间正在烦恼做计算所花费的时间太长，正在徘徊于是否要去进行CUDA编程的学习，这时忽然我关注到了Taichi语言，此处需要非常感谢我的好哥们丁一，正是他带我关注到了他，虽然Taichi仍然很年轻，但是由于他简单易学的语法，使得我在编程方式的转换上十分轻松，不需要去学习大量的计算机基础知识，也不需要去学习复杂的CUDA编程，只需要在简单适应并行编程的思维逻辑之后，用ti.init(arch=ti.gpu)，编译运行之后就可以实现高效的并行计算，而不用去苦等好几个小时只为一个模拟结果；由于Taichi自身的开源的属性，使得我在学习他的过程中可以省去很多后顾之忧，并且可以和开发者直接提出自己所遇到的问题，并且可以提出自己的想法，头脑风暴的过程就是自己提高的过程，每次与数值计算群里的伙伴们的交流，都使得我获益良多；Taichi自身是嵌入在Python中的一个库，因此可以无缝与python的其他各个库联合使用，使得我可以在计算的过程中实时监测流场中的数据变化，极大的增强了我的数据处理能力，并且为我闭环优化所研究的系统提供了可能，对我的科研工作帮助极大。综上所述，这就是我将LBM与Taichi强强联合的原因，我也会尽力为Taichi语言的发展和Taichi社区的繁荣奉献自己的一份力量。

—— \tau 于2023.1.30记

[Tin]

写的真好。Taichi做着"flattern the learning curve"的工作，而这篇文章，整个社区也在做着同样的事情。真好。