从格子Boltzmann方程到宏观流体问题的常用还原方法

这是我之前的一次期末大报告作业，里面的内容也算是从其他文章里搬运的，在此重新誊抄出来。内容如有错误，欢迎各位指正。
[注] 下文提到的 $x$ 、 $u$ 和 $\xi$ 无论有没有加粗均为向量。

本文同样也发在了知乎专栏上。

一、 单松弛格子Boltzmann方程

1.1 Boltzmann方程

格子Boltzmann方法基于Boltzmann方程（即(1)式），该方程引入速度分布函数，描述了非热力学平衡状态的热力学系统统计行为：

$$\frac{\partial f}{\partial t}+\xi \cdot \frac{\partial f}{\partial x}+\frac{F}{m}\cdot \frac{\partial f}{\partial \xi }=\Omega (f)\text{\hspace{1em}(1)}$$

(1)式中$\Omega (f)$为碰撞项， $\xi =\frac{\partial x}{\partial t}$ 。速度分布函数的定义是：在时刻$t$，以$x$为中心的空间微元 $\mathrm{d}x$ 内，速度在 $\xi$ 和 $(\xi +\mathrm{d}\xi )$ 之间的分子的数量为 $f\mathrm{d}x\mathrm{d}\xi$ 。记分子质量为 $m$ ，流体宏观速度为$u$，单位体积内能为$e$，分子数密度为$n$，则：

$$\rho =mn=m\int f\mathrm{d}\xi$$

$$\rho u=m\int \xi f\mathrm{d}\xi$$

$$\rho E=\rho e+\frac{1}{2}\rho u^2=m\int f\frac{{\xi }^2}{2}\mathrm{d}\xi$$

碰撞项 描述了分子碰撞对整个系统运动产生的影响，它具有两个特点：①满足质量、动量和能量守恒；②能够反映系统趋于平衡态的趋势。Bhatnagur、Gross和Krook三人在此基础上提出了最为简单且著名的BGK模型，即(2)式。

$$\Omega (f)=-\frac{1}{{\tau }_c}\left(f-f^{(eq)}\right)\text{\hspace{1em}(2)}$$

其中${\tau }_c$为松弛时间， $f^{(eq)}$ 为平衡态分布函数。 可根据Boltzmann H定理计算得：

$$f^{(eq)}=\frac{n}{(2\pi RT{)}^{3/2}}\exp \left(-\frac{(\xi -\mathbf{u}{)}^2}{2RT}\right)\text{\hspace{1em}(3)}$$

1.2 格子Boltzmann方程

为了将Boltzmann方程应用于实际的数值模拟，需要将流体视为大量的离散粒子。离散粒子于流体分子不同，它们驻留在网格上，并按一定的规则在网格上发生碰撞和迁移。因此，1.2节将 $mf$ 记为 $f$，此时 $f$ 为密度分布函数。对密度分布函数 $f$ ，不考虑 $F$ 的Boltzmann-BGK方程为式。

$$\frac{\partial f}{\partial t}+\xi \cdot \frac{\partial f}{\partial x}=-\frac{1}{{\tau }_c}\left[f-f^{(eq)}\right]\text{\hspace{1em}(4)}$$

其宏观量按如下方式确定

$$\rho =\int f\mathrm{d}\xi ，\rho \mathrm{u}=\int \xi \mathrm{f}\mathrm{d}\xi ，\rho \mathrm{E}=\int \mathrm{f}\frac{{\xi }^2}{2}\mathrm{d}\xi$$

对 $f$ 和其对应的 $f^{(eq)}$ 可进行如下离散：$f_i=w_{i}f$ ， $f_i^{(eq)}=w_i{f}^{(eq)}$ 。其中下标$i$表示特征方向${\mathbf{c}}_i$， $f_i$ 和$f_i^{(eq)}$均沿着${\mathbf{c}}_i$方向运动。据此，沿着特征方向离散(4)式得到的差分格式为

$$f_i(x+{\xi }_i\Delta t,t+\Delta t)-f_i(x,t)=-\frac{1}{\tau }\left[f_i(x,t)-f_i^{(eq)}(x,t)\right]\text{\hspace{1em}(5)}$$

其中 $\tau =\frac{{\tau }_c}{\Delta t}$ 为无量纲松弛时间。由于对应的$f^{(eq)}$为：

$$f^{(eq)}=\rho \cdot \exp (-\frac{(\xi -u{)}^2}{2RT})\text{\hspace{1em}(6-1)}$$

所以将$f^{(eq)}$展开至2阶并离散，可得：

$$f_i^{(eq)}=w_i\rho \left[1+\frac{{\xi }_i\cdot u}{c_s^2}+\frac{({\xi }_i\cdot u{)}^2}{2c_s^4}-\frac{u\cdot u}{2c_s^2}\right]\text{\hspace{1em}(6-2)}$$

其中 $c_s=\sqrt{RT}$ 在LBM模拟中一般被称作格子声速。

二、Chapman-Enskog展开分析

将基本碰撞不变量 ${\phi }_i=1,\xi ,\frac{{\xi }^2}{2}$ 乘以(1)式两端，并对 $\xi$ 积分有：

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

$$\frac{\partial (\rho u)}{\partial t}+\nabla \cdot (\rho \mathbf{uu}+P)=\rho \mathbf{a}\text{\hspace{1em}(7-2)}$$

$$\frac{\partial (\rho E)}{\partial t}+\nabla \cdot (\rho \mathbf{u}E+P\cdot \mathbf{u}+\mathbf{Q})=\rho \mathbf{a}\cdot \mathbf{u}\text{\hspace{1em}(7-3)}$$

Chapman-Enskog展开（下文简称C-E展开）是处理Boltzmann方程的一种分析方法，该方法引入小量$\epsilon$作为展开因子（一般与Knudsen数同阶），将BGK方程在不同的尺度上展开，即：

$$f=\sum _{n=0}^{\infty }{\epsilon }^n{f}^{(n)}\text{\hspace{1em}(8)}$$

并将分布函数、导数、物理量等都按照$\epsilon$的不同阶次展开。该方法对 $f^{(n)}$ 进行了限制，使其满足：$\int f^{(n)}{\phi }_i\mathrm{d}\xi =0$， ${\phi }_i(i=1,..,5)$ 为基本碰撞不变量。将 $f$ 视为 $\rho$ 、$\mathbf{u}$ 和 $T$ 及其各阶导数的泛函，而时间偏导项则展开为

$$\frac{\partial }{\partial t}=\sum _{n=0}^{\infty }\frac{{\partial }_n}{\partial t}\text{\hspace{1em}(9)}$$

将(8)式和(9)式代入到(7)式中。得 $O({\epsilon }^0)$ 方程组为：

$$\frac{{\partial }_0\rho }{\partial t}+\nabla \cdot (\rho \mathbf{u})=0\text{\hspace{1em}(10-1)}$$

$$\frac{{\partial }_0(\rho u)}{\partial t}+\nabla \cdot (\rho \mathbf{uu}+P^{(0)})=\rho \mathbf{a}\text{\hspace{1em}(10-2)}$$

$$\frac{{\partial }_0(\rho E)}{\partial t}+\nabla \cdot (\rho \mathbf{u}E+P^{(0)}\cdot \mathbf{u}+{\mathbf{Q}}^{(0)})=\rho \mathbf{a}\cdot \mathbf{u}\text{\hspace{1em}(10-3)}$$

以及 $O({\epsilon }^n)$ 方程组为 $(n>0)$ ：

$$\frac{{\partial }_n\rho }{\partial t}=0\text{\hspace{1em}(11-1)}$$

$$\frac{{\partial }_n(\rho u)}{\partial t}+\nabla \cdot P^{(n)}=0\text{\hspace{1em}(11-2)}$$

$$\frac{{\partial }_n(\rho E)}{\partial t}+\nabla \cdot (P^{(n)}\cdot \mathbf{u}+{\mathbf{Q}}^{(n)})=0\text{\hspace{1em}(11-3)}$$

其中
$$C=\mathbf{\xi }-\mathbf{u},P^{(n)}=m\int CC{f}^{(n)}\mathrm{d}\xi ,Q^{(n)}=m\int C\frac{C^2}{2}{f}^{(n)}\mathrm{d}\xi$$

只需取精确到$O({\epsilon }^1)$的宏观方程组进行分析，即可还原出Naiver-Stokes方程。且具有如下关系：

$$P=P^{(0)}+\epsilon P^{(1)}=p\mathbf{I}-2\epsilon \mu (\mathbf{S}-\frac{\mathrm{Tr}(S)}{3}\mathbf{I})$$

$$Q=Q^{(0)}+\epsilon Q^{(1)}=-\epsilon \kappa \nabla T$$

其中：$\kappa$为热传导系数， $\mu$ 为粘性系数。

离散的LBGK方程也可采用类似的方法展开，还原出宏观流体运动方程（见附录A）。展开方式如下：

$$f_i=f_i^{(0)}+\epsilon f_i^{(1)}+{\epsilon }^2{f}_i^{(2)}+...$$

$$\frac{\partial }{\partial x}=\epsilon \frac{\partial }{\partial x_1},x_1=\epsilon x,\frac{\partial }{\partial t}=\epsilon \frac{\partial }{\partial t_1}+{\epsilon }^2\frac{\partial }{\partial t_2},F_i=\epsilon F_i^{(1)}$$

略去$O({\epsilon }^3)$项，并对$O({\epsilon }^0)$项、$O({\epsilon }^1)$项和$O({\epsilon }^2)$项求0阶和1阶速度矩，即可得到不可压缩流体的Navier-Stokes方程。由于在分析过程中略去了$O({\epsilon }^3)$项和速度的3阶项，导致结果缺失了部分物理信息。

三、 Grad-13矩方法及其发展

3.1 Grad的矩分析方法

Grad -13矩是一种由数学家Harold Grad提出的分析方法，其基于Hermite多项式进行分析。对于一个长度为 $N$ 的向量 $\mathbf{x}$ ，若定义权函数

$$\omega (\mathbf{x})=\frac{1}{(2\pi {)}^{N/2}}\exp \left(-\frac{\mathbf{x}\cdot \mathbf{x}}{2}\right)\text{\hspace{1em}(12)}$$

则其对应的$n$阶Hermite多项式写作

$${\mathcal{H}}^{(n)}(\mathbf{x})=\frac{(-1{)}^n}{\omega (\mathbf{x})}{\nabla }^n\omega (\mathbf{x})\text{\hspace{1em}(13)}$$
由Hermite多项式的正交性，可证明得到：

$${\int }_{-\infty }^{+\infty }\omega (\mathbf{x})H_i^{(m)}(\mathbf{x})H_j^{(n)}(\mathbf{x})\mathrm{d}\mathbf{x}={\delta }_{\mathrm{mn}}{\delta }_{\mathrm{ij}}\text{\hspace{1em}(14)}$$

据此，可将任意一个函数 $f(x)$ 使用Hermite多项式进行展开¹。即：

$$f(\mathbf{x})={\omega }^{\frac{1}{2}}\sum _{n=0}^{\infty }\sum _{i=1}^N{a}_i^{(n)}{H}_i^{(n)}(\mathbf{x})\text{\hspace{1em}(15)}$$
所以

$${\int }_{-\infty }^{-\infty }{\omega }^{1/2}f(x){\mathcal{H}}_j^{(m)}(x)\mathrm{d}x=\sum _{i=1}^N{a}_i^{(m)}{\delta }_{ij}^m=m!a_i^m\text{\hspace{1em}(16)}$$

其中 $a_i^{(m)}$ 可借由Hermite多项式的正交性进行求解。 ${\delta }_{ij}^m$ 是一个记号，表示 $m$个 ${\delta }_{ij}$的乘积之和， $i$ 和 $j$ 分别来自下标集 $(i_1,i_2,...,i_m)$ 和 $(j_1,j_2,...,j_m)$ 。
在此数学基础之上，Harold Grad将分布函数 $f$ 展开为²：

$$f(x,\xi ,t)=\omega \sum _{n=0}^{\infty }\frac{1}{n!}{a}^{(n)}{H}^{(n)}(\xi )\text{\hspace{1em}(17)}$$

根据式(16)， $a^{(n)}$的表达式为：

$$a^{(n)}={\int }_{-\infty }^{+\infty }f(x,\xi ,t)H^{(n)}(\xi )d\xi \text{\hspace{1em}(18)}$$

可得$a^{(n)}$的前4项为：

$$a^{(0)}={\int }_{-\infty }^{+\infty }f(x,\xi ,t)d\xi =\rho \text{\hspace{1em}(19-1)}$$

$$a^{(1)}={\int }_{-\infty }^{+\infty }f(x,\xi ,t)\xi d\xi =\rho u\text{\hspace{1em}(19-2)}$$

$$a^{(2)}={\int }_{-\infty }^{+\infty }f(x,\xi ,t)({\xi }^2-\delta )d\xi =\rho ({\mathbf{u}}^2-\delta )+\mathbf{P}\text{\hspace{1em}(19-3)}$$

$$a^{(3)}={\int }_{-\infty }^{+\infty }f(x,\xi ,t)({\xi }^3-\xi \delta )d\xi =(D-1)\rho {\mathbf{u}}^3+\mathbf{u}{a}^{(2)}+\mathbf{Q}\text{\hspace{1em}(19-4)}$$

注意到从$a^{(0)}$到 $a^{(3)}$ 都描述了流场的宏观量，所以通过上述表达式的值可以解出流场宏观量。与C-E分析方法不同的是分析过程中没有舍去高阶量。因此这种展开方法得到的式子合理地引入了更高阶的矩，物理意义更加明确。

3.2 Gauss-Hermite积分的引入和离散化BGK方程

Grad的思路开创了描述的新方式，而单肖文、袁学锋、陈沪东三位学者在此基础上进一步发展，于2006年在 Journal of Fluid Mechanics 上提出了将原思路拓展到离散的Boltzmann-BGK方程中的方法³，简化了求解过程。
首先，基于式(16)将 $f$ 截断至前 $N$ 阶，得：

$$f\approx f^N=\omega \sum _{n=0}^N\frac{1}{n!}{a}^{(n)}{H}^{(n)}(\xi )\text{\hspace{1em}(20)}$$

$f^N$ 和 $f$ 有着完全相同的前 $N$ 阶矩。通过截断操作将BGK方程对宏观流动方程的近似保持在动力学层面而不是热力学层面。此时$f^N$可表示为：

$$f^N(x,\xi ,t)H^{(n)}(\xi )=\omega (\xi )p(x,\xi ,t)\text{\hspace{1em}(21)}$$

其中 $p(x,\xi ,t)$ 为阶数不超过 $2N$ 的多项式。使用Gauss-Hermite积分可将 $a^{(n)}$ 精确地展开为 $p(x,\xi ,t)$ 的加权和。

$$a^{(n)}=\int \omega (\xi )p(x,\xi ,t)d\xi =\sum _{a=1}^d{w}_{a}p(x,{\xi }_a,t)=\sum _{a=1}^d\frac{w_a}{\omega ({\xi }_a)}{f}^N(x,{\xi }_a,t)H^{(n)}({\xi }_a)\text{\hspace{1em}(22)}$$

式(22)中 $w_a$ 和 ${\mathbf{\xi }}_a$ 是Gauss-Hermite积分的权重和特征向量（ $a=1,2,...,d$ ）。所以使用 $f^N(x,{\mathbf{\xi }}_a,t)$ 表示 $f^N$ 的前 $N$ 阶矩在数学上可行。
据此可定义离散化的分布函数 $f_a$ （ $a=1,2,...,d$ ），其对应的权重和特征方向为 $w_a$ 和 ${\xi }_a$ ：

$$f_a(x,t)=\frac{w_a}{\omega ({\mathbf{\xi }}_a)}f(x,{\mathbf{\xi }}_a,t)\text{\hspace{1em}(23)}$$

则宏观量可表达为：

$$\rho =\sum _{a=1}^d{f}_a,\rho \mathbf{u}=\sum _{a=1}^d{f}_a{\xi }_a,\rho \mathbf{uu}+\mathbf{P}=\sum _{a=1}^d{f}_a{\xi }_a{\xi }_a\text{\hspace{1em}(24)}$$

下面将(1)式的碰撞项换为(2)式，并把 $\frac{\mathbf{F}}{m}\frac{\partial f}{\partial \xi }$ 直接记作外力项 $-\mathbf{F}(\xi )$ ，得

$$\frac{\partial f}{\partial t}+\xi \cdot \frac{\partial f}{\partial x}=\frac{-1}{\tau }(f-f^{(eq)})+\mathbf{F}(\xi )\text{\hspace{1em}(25)}$$

对式(25)，令 $\mathbf{\xi }={\mathbf{\xi }}_a$ ，并在两端乘上 $ \frac{w_a}{\omega( \boldsymbol{\xi}_a)}$ ，即可得到使用 $f_a$ 离散的方程（ $a=1,2,...,d$ ）：

$$\frac{\partial f_a}{\partial t}+{\xi }_a\cdot \nabla f_a=-\frac{1}{\tau }(f_a-f_a^{(eq)})+F_a\text{\hspace{1em}(26-1)}$$

$$f_a^{(eq)}=\frac{w_a}{\omega ({\xi }_a)}{f}^{(eq)}({\xi }_a)\text{\hspace{1em}(26-2)}$$

$$F_a=\frac{w_a}{\omega ({\xi }_a)}\mathbf{F}({\xi }_a)\text{\hspace{1em}(26-3)}$$

此时只需将 $f^{(eq)}$ 和 $F$ 进行类似的展开，即可完成对整个系统的离散。

附录A． 使用C-E展开方法将离散的LBGK方程还原为Navier-Stokes方程

LBE方程为(A-1.1)式。源项 $\mathbf{F}$ 和 $f^{(\mathrm{eq})}$ 为由郭照立、郑楚光和施保昌提出的格式，即(A-1.2)式和(A-1.3)式。[注意，这里的 $f$ 是密度分布函数]

$$f_i(x+{\xi }_i\Delta t,t+\Delta t)-f_i(x,t)=-\frac{1}{\tau }\left[f_i(x,t)-f_i^{(eq)}(x,t)\right]+\Delta t\cdot F_i(x,t)\text{\hspace{1em}(A-1.1)}$$

$$F_i=(1-\frac{\Delta t}{2\tau })w_i\left[\frac{{\xi }_{\mathbf{i}}-\mathbf{u}}{c_s^2}+\frac{({\xi }_{\mathbf{i}}\cdot \mathbf{u}){\xi }_{\mathbf{i}}}{c_s^4}\right]\cdot F\text{\hspace{1em}(A-1.2)}$$

$$f_i^{(eq)}(x,t)=\rho w_i\left[1+\frac{{\xi }_{\mathbf{i}}\cdot \mathbf{u}}{c_s^2}+\frac{\mathbf{uu}:({\xi }_{\mathbf{i}}{\xi }_{\mathbf{i}}-c_s^2\mathbf{I})}{2c_s^4}\right]\text{\hspace{1em}(A-1.3)}$$

并且对于 $f_i$ 和 $f_i^{(eq)}$ 有：

$$\{\begin{array}{l}\rho ={\sum }_i{f}_i\\ \rho \mathbf{u}={\sum }_i{\xi }_{\mathbf{i}}{f}_i+\frac{\Delta t}{2}F\end{array}$$

和

$$\{\begin{array}{l}\rho ={\sum }_i{f}_i^{(eq)}\\ \rho \mathbf{u}={\sum }_i{\xi }_{\mathbf{i}}{f}_i^{(eq)}\\ \rho (\mathbf{uu}+c_s^2\mathbf{I})={\sum }_i{\xi }_{\mathbf{i}}{\xi }_{\mathbf{i}}{f}_i^{(eq)}\end{array}$$

将(A-1.1)式在 $(\mathbf{x},t)$处展开，取至 2 阶项得：
$$\begin{array}{rl}\Delta t\frac{\partial f_i}{\partial t}+{\xi }_i\Delta t\cdot \frac{\partial f_i}{\partial x}& +\frac{1}{2}\left[\Delta t^2\frac{{\partial }^2{f}_i}{\partial t^2}+2\Delta t({\xi }_i\Delta t)\cdot \frac{{\partial }^2{f}_i}{\partial t\partial x}+\Delta t^2({\xi }_i{\xi }_i):\frac{{\partial }^2{f}_i}{\partial x^2}\right]\\ & +\frac{1}{\tau }\left[f_i-f_i^{(eq)}\right]-\Delta t{F}_i=0\end{array}\text{\hspace{1em}(A-2)}$$

代入多尺度展开等式，并舍去 ${\epsilon }^3$ 项和更高阶的项，可得：
$$\begin{array}{rl}(\epsilon \frac{\partial f_i^{(0)}}{\partial t_1}+{\epsilon }^2\frac{\partial f_i^{(1)}}{\partial t_1}& +{\epsilon }^2\frac{\partial f_i^{(0)}}{\partial t_2})+{\xi }_i\cdot (\epsilon \frac{\partial f_i^{(0)}}{\partial x_1}+{\epsilon }^2\frac{\partial f_i^{(1)}}{\partial x_1})\\ & +\frac{\Delta t}{2}\left[{\epsilon }^2\frac{{\partial }^2{f}_i^{(0)}}{\partial t_1^2}+2{\xi }_i{\epsilon }^2\cdot \frac{{\partial }^2{f}_i^{(0)}}{\partial t_1\partial x_1}+({\xi }_i{\xi }_i):({\epsilon }^2\frac{{\partial }^2{f}_i^{(0)}}{\partial x^2})\right]\\ & +\frac{1}{\tau \Delta t}\left[f_i^{(0)}+\epsilon f_i^{(1)}+{\epsilon }^2{f}_i^{(2)}-f_i^{(eq)}\right]-\epsilon F_i^{(1)}+O({\epsilon }^3)=0\end{array}\text{\hspace{1em}(A-3)}$$

从(A-3)式中提取$O({\epsilon }^0)$项、$O({\epsilon }^1)$项和$O({\epsilon }^2)$项，并令 $D_k=\frac{\partial }{\partial t_k}+{\xi }_k\cdot \frac{\partial }{\partial {\mathbf{x}}_k}$ 。(A-3)式对任意 $\epsilon$ 都成立，则关于各阶的系数均为0，即：
$$O({\epsilon }^0)=\frac{1}{\tau \Delta t}(f_i^{(0)}-f_i^{(eq)})=0\text{\hspace{1em}(A-4.1)}$$
$$O({\epsilon }^1)=D_1{f}_i^{(0)}+\frac{f_i^{(1)}}{\tau \Delta t}-F_i^{(1)}=0\text{\hspace{1em}(A-4.2)}$$
$$O({\epsilon }^2)=\frac{\partial f_i^{(0)}}{\partial t_2}+(1-\frac{1}{2\tau })D_1{f}_i^{(1)}+\frac{f_i^{(2)}}{\tau \Delta t}+\frac{\Delta t}{2}{D}_1{F}_i^{(1)}=0\text{\hspace{1em}(A-4.3)}$$

根据多尺度展开的变量关系，可得：
${\sum }_i{f}_i^{(1)}={\sum }_i{f}_i^{(2)}=0$ ， ${\sum }_i{\xi }_i{f}_i^{(2)}=0$ ， ${\sum }_i{\xi }_i{f}_i^{(1)}=-\frac{\Delta t}{2}{F}_i^{(1)}$

$$\{\begin{array}{c}{\sum }_i{F}_i=0\\ {\sum }_i{\xi }_i{F}_i=(1-\frac{1}{2\tau })F\\ {\sum }_i{\xi }_i{\xi }_i{F}_i=(1-\frac{1}{2\tau })2uF\end{array}$$
$$\{\begin{array}{c}{\sum }_i{F}_i^{(1)}=0\\ {\sum }_i{\xi }_i{F}_i^{(1)}=(1-\frac{1}{2\tau })F^{(1)}\\ {\sum }_i{\xi }_i{\xi }_i{F}_i^{(1)}=(1-\frac{1}{2\tau })2u{F}^{(1)}\end{array}$$

为方便下文分析，记： ${\Pi }^{(k)}={\sum }_i{\xi }_i{\xi }_i{f}_i^{(k)}$ ， ${\Gamma }^{(k)}={\sum }_i{\xi }_i{\xi }_i{\xi }_i{f}_i^{(k)}$ 。

对$O({\epsilon }^1)$项求1阶和2阶速度矩，得：

$$\frac{\partial \rho }{\partial t_1}+\frac{\partial (\rho u)}{\partial x_1}=0\text{\hspace{1em}(A-5.1)}$$

$$\frac{\partial (\rho u)}{\partial t_1}+\frac{\partial {\Pi }^{(0)}}{\partial x_1}=0\text{\hspace{1em}(A-5.2)}$$

$${\Pi }^{(0)}=\rho (\mathbf{uu}+{\mathrm{c}}_{\mathrm{s}}^2\mathbf{I})\text{\hspace{1em}(A-5.3)}$$

对$O({\epsilon }^2)$项求1阶和2阶速度矩，得：

$$\frac{\partial \rho }{\partial t_2}=0\text{\hspace{1em}(A-6.1)}$$

$$\frac{\partial (\rho \mathbf{u})}{\partial t_2}+(1-\frac{1}{2\tau })\frac{\partial {\Pi }^{(1)}}{\partial {\mathbf{x}}_1}+\Delta t(1-\frac{1}{2\tau })\frac{\partial {\mathbf{F}}^{(1)}}{\partial {\mathbf{x}}_1}=0\text{\hspace{1em}(A-6.2)}$$

由(A-4.2)式，得： f i ( 1 ) = τ Δ t { D 1 f i ( 0 ) − F i ( 1 ) } f_i^{(1)} = \tau \Delta t \{ D_1 f_i^{(0)} - F_i^{(1)}\} fi(1)​=τΔt{D1​fi(0)​−Fi(1)​} ，所以

$${\Pi }^{(1)}=\sum _i{\mathbf{\xi }}_i{\mathbf{\xi }}_i{f}_i^{(1)}=-\tau \Delta t[\frac{\partial (\rho \mathbf{u}\mathbf{u})}{\partial t_1}+c_s^2\frac{\partial (\rho \mathbf{I})}{\partial t_1}+\frac{\partial {\Gamma }^{(0)}}{\partial {\mathbf{x}}_1}-(1-\frac{1}{2\tau })2\mathbf{u}{\mathbf{F}}^{(1)}]\text{\hspace{1em}(A-7)}$$

$\mathbf{u}{\mathbf{F}}^{(1)}$是Grad的论文中的记号， $({uF}^{(1)}{)}_{lm}=\frac{1}{2}(u_l{F}_m^{(1)}+u_m{F}_l^{(1)})$ 。下文为方便分析，取 ${\Pi }_{lm}^{(1)}$ 分析。

对第1项，有： $\frac{\partial (\rho u_l{u}_m)}{\partial t_1}=\rho u_l\frac{\partial u_m}{\partial t_1}+u_m\frac{\partial (\rho u_l)}{\partial t_1}$ 。由(A-5.2)式，可得⁴：

$$\rho \frac{\partial u_m}{\partial t_1}+u_m\frac{\partial \rho }{\partial t_1}+\frac{\partial }{\partial x_{1n}}\left(\rho u_m{u}_n+\rho c_s^2{\delta }_{mn}\right)=F_m^{(1)}\text{\hspace{1em}(A-8.1)}$$

$$\frac{\partial (\rho u_l)}{\partial t_1}+\frac{\partial }{\partial x_{1n}}\left(\rho u_l{u}_n+\rho c_s^2{\delta }_{ln}\right)=F_l^{(1)}\text{\hspace{1em}(A-8.2)}$$

并且

$$\frac{\partial (\rho u_l{u}_m{u}_n)}{\partial x_{1n}}=-u_l{u}_m\frac{\partial (\rho u_n)}{\partial x_{1n}}+u_l\frac{\partial (\rho u_m{u}_n)}{\partial x_{1n}}+u_m\frac{\partial (\rho u_l{u}_n)}{\partial x_{1n}}\text{\hspace{1em}(A-8.3)}$$

所以
$$\frac{\partial (\rho u_l{u}_m)}{\partial x_{1n}}=(u_l{F}_m^{(1)}+u_m{F}_l^{(1)})-\frac{\partial (\rho u_l{u}_m{u}_n)}{\partial x_{1n}}-c_s^2\left(u_l\frac{\partial \rho }{\partial x_{1m}}+u_m\frac{\partial \rho }{\partial x_{1l}}\right)\text{\hspace{1em}(A-8.4)}$$

所以第2项的分量为： $c_s^2{\delta }_{lm}\frac{\partial \rho }{\partial t_1}$ 。

对第3项的分量，由于
$${\Gamma }_{lmn}^{(0)}=\rho c_s^2(u_l{\delta }_{mn}+u_n{\delta }_{lm}+u_m{\delta }_{ln})$$
，所以：

$$\frac{\partial {\Gamma }_{lmn}^{(0)}}{\partial x_{1n}}=c_s^2({\delta }_{mn}\frac{\partial \rho u_l}{\partial x_{1n}}+{\delta }_{lm}\frac{\partial \rho u_n}{\partial x_{1n}}+{\delta }_{ln}\frac{\partial \rho u_m}{\partial x_{1n}})\text{\hspace{1em}(A-8.5)}$$

所以

$${\Pi }_{lm}^{(1)}=-\tau \Delta t\left[-\frac{\partial (\rho u_l{u}_m{u}_n)}{\partial x_{1n}}+\frac{1}{2\tau }(u_l{F}_m^{(1)}+u_m{F}_l^{(1)})+\rho c_s^2(\frac{\partial u_l}{\partial x_{1m}}+\frac{\partial u_m}{\partial x_{1l}})\right]\text{\hspace{1em}(A-8.6)}$$

略去${\Pi }_{lm}^{(1)}$ 中的速度的3阶项（即 ∂ ∂ x 1 n { ρ u l u m u n } \frac{\partial}{\partial x_{1n}} \{ \rho u_l u_m u_n \} ∂x1n​∂​{ρul​um​un​} ），有

$${\Pi }^{(1)}=-\Delta t\cdot {\mathbf{uF}}^{(1)}-\tau \rho c_s^2\Delta t\left[\frac{\partial \mathbf{u}}{\partial {\mathbf{x}}_1}+{\left(\frac{\partial \mathbf{u}}{\partial {\mathbf{x}}_1}\right)}^{\mathrm{T}}\right]\text{\hspace{1em}(A-9)}$$
回代（A-6.2）式，得等价形式为

$$\frac{\partial (\rho \mathbf{u})}{\partial t_2}+c_s^2\Delta t(\tau -\frac{1}{2})\frac{\partial }{\partial {\mathbf{x}}_1}\left\{\rho \left[\frac{\partial \mathbf{u}}{\partial {\mathbf{x}}_1}+{\left(\frac{\partial \mathbf{u}}{\partial {\mathbf{x}}_1}\right)}^{\mathrm{T}}\right]\right\}=0\text{\hspace{1em}(A-10)}$$

将(A-10)式和(A-5.2)式联立，得

∂ ( ρ u ) ∂ t + ∂ ( ρ u u ) ∂ x = − ∂ ( ρ c s 2 ) ∂ x + c s 2 Δ t ( τ − 1 2 ) ∂ ∂ x { ρ [ ∂ u ∂ x + ( ∂ u ∂ x ) T ] } + F (A-11) \frac{\partial (\rho u)}{\partial t} + \frac{\partial (\rho u u)}{\partial x} = - \frac{\partial (\rho c_s^2)}{\partial x} + c_s^2 \Delta t (\tau - \frac{1}{2}) \frac{\partial}{\partial x}\{\rho [\frac{\partial u}{\partial x} + (\frac{\partial u}{\partial x})^T ]\} + F \tag{A-11} ∂t∂(ρu)​+∂x∂(ρuu)​=−∂x∂(ρcs2​)​+cs2​Δt(τ−21​)∂x∂​{ρ[∂x∂u​+(∂x∂u​)T]}+F(A-11)

此时令 $p=c_s^2\rho$ ， $\nu =c_s^2\Delta t(\tau -\frac{1}{2})$ ，即可还原回Navier-Stokes方程

附录B．Gauss-Hermite积分

对于任意一维函数 $f(\xi )$ ，高斯积分通过寻找积分点 ${\xi }_k$ 和其对应权重 $w_k$ 实现对 ${\int }_a^b\omega (\xi )f(\xi )d\xi$ 的近似，即

$${\int }_a^b\omega (\xi )f(\xi )d\xi \approx \sum _{k=1}^n{w}_k{\xi }_k\text{\hspace{1em}(B-1)}$$

其中 $\omega (\xi )$ 是任意的权重函数。n 点高斯求积的积分点正是对应的第n个正交多项式的根，此时 ${\xi }_k$ 对应的权重 $w_k$ 为

$$w_k=\frac{<P_{n-1},P_{n-1}>}{P_{n-1}({\xi }_k)P_n^{{}^{\prime }}({\xi }_k)}\text{\hspace{1em}(B-2)}$$

其中 $P_n^{{}^{\prime }}=\frac{d{P}_n}{dx}$ 。(B-2)式的精度为 $2n-1$ 。
所以在一维Gauss-Hermite积分中，取 $P_n={\mathcal{H}}^n$ ， ${\mathcal{H}}^n$ 的权函数和表达式见(12)式和(13)式（下文同理）。 n点高斯求积的积分点即为 ${\mathcal{H}}^n$ 的根。因为 $\frac{d{\mathcal{H}}^n}{d\xi }=\xi {\mathcal{H}}^n-{\mathcal{H}}^{n+1}=n{\mathcal{H}}^{n-1}$ ，所以有

$$w_k=n!/[n{\mathcal{H}}^{(n-1)}({\xi }_k){]}^2\text{\hspace{1em}(B-3)}$$

对于更高维的情况，可做类似构造³。分析如下积分

$$\frac{1}{(2\pi {)}^{D/2}}\int \exp (-\frac{\xi \cdot \xi }{2})p(\xi )\mathrm{d}\xi \text{\hspace{1em}(B-4)}$$

其中 $\xi$ 是长度为D的向量。 $p(\xi )$ 是D维n次多项式，可写为(B-5)式

$$p(\xi )=\sum _{n_1+n_2+...+n_D\le n}{a}_{n_1{n}_2...n_D}\prod _{j=1}^D{\xi }_j^{n_j}\text{\hspace{1em}(B-5)}$$

其中 $a_{n_1{n}_2...n_D}$为常数。

考虑(B-5)式中的任意一项 $\prod _{j=1}^D{\xi }_j^{n_j}$ （由于$a_{n_1{n}_2...n_D}$ 为常数所以可暂不处理，在最后一步被归并到常数项中），回代到式(B-4)中，得

$$\frac{1}{(2\pi {)}^{D/2}}\int \exp (-\frac{\xi \cdot \xi }{2})p(\xi )\mathrm{d}\xi =\prod _{j=1}^D\left(\frac{1}{\sqrt{2\pi }}\int \exp (-\frac{{\xi }_j^2}{2}){\xi }_j^{n_j}\mathrm{d}{\xi }_j\right)=\prod _{j=1}^D(\sum _{k=1}^n{w}_k{\xi }_k)\text{\hspace{1em}(B-6.1)}$$

其中 $w_k$和 ${\xi }_k$是一维n次多项式高斯积分的权重和积分点。积分结果被转化为多个一维Gauss-Hermite积分的连乘。注意到(B-6.1)式等价于

$$\frac{1}{(2\pi {)}^{D/2}}\int \exp (-\frac{\xi \cdot \xi }{2})p(\xi )\mathrm{d}\xi =\sum _{k_1=1}^n...\sum _{k_D=1}^n[(w_{k_1}{w}_{k_2}...w_{k_D})({\xi }_{k_1}^{n_1}{\xi }_{k_2}^{n_2}...{\xi }_{k_D}^{n_D})]\text{\hspace{1em}(B-6.2)}$$

如果定义 ${\xi }_{k_1...k_D}=({\xi }_{k_1},{\xi }_{k_2},...,{\xi }_{k_D})$ 和 $w_{k_1...k_D}=w_{k_1}{w}_{k_2}...w_{k_D}$ ，那么(B-4)式等价为

$$\frac{1}{(2\pi {)}^{D/2}}\int \exp (-\frac{\xi \cdot \xi }{2})p(\xi )\mathrm{d}\xi =\sum w_{k_1...k_D}p({\xi }_{k_1...k_D})\text{\hspace{1em}(B-7)}$$

此时则将(B-4)式转换为与一维Gauss-Hermite积分类似的形式，即(B-7)式。

参考

---

1. GRAD H. Note on N-Dimensional Hermite Polynomials[J]. Communications on Pure and Applied Mathematics, 1949, 2(4): 325–330. http://doi.org/10.1002/cpa.3160020402 ↩︎

2. GRAD H. On the Kinetic Theory of Rarefied Gases[J]. Communications on Pure and Applied Mathematics, 1949, 2(4): 331–407. http://doi.org/10.1002/cpa.3160020403 ↩︎

3. SHAN X, YUAN X-F, CHEN H. Kinetic Theory Representation of Hydrodynamics: A Way beyond the Navier–Stokes Equation[J]. Journal of Fluid Mechanics, 2006, 550(1): 413. http://doi.org/10.1017/S0022112005008153 ↩︎ ↩︎

4. CAO W. Investigation of the Applicability of the Lattice Boltzmann Method to Free-Surface Hydrodynamic Problems in Marine Engineering[D/OL]. Laboratoire de recherche en Hydrodynamique, Énergétique et Environnement Atmosphérique (LHEEA): École centrale de Nantes, 2019. https://tel.archives-ouvertes.fr/tel-02383174 ↩︎