0 引　言

诞生于格子气自动机^[1]的格子Boltzmann方法（LBM），与传统直接求解宏观量的计算流体力学（CFD）方法相比，由于避开了Navier-Stokes方程非线性项的处理，且没有求解压力泊松方程这一繁重过程，使得LBM具有程序编制简单、并行计算效率高的优点^[2]，在湍流、多相流、多孔介质流动及渗流、颗粒流、微流动、气动声学等领域的模拟中得到了广泛应用^[3-6]。早期的LBM离散模型都是基于半经验理论推导，其约束条件往往局限于满足质量和动量守恒方程，并以低马赫数假设为基础，能量方程无法正确还原，使得LBM方法只适用于等温低速弱可压流动成为了较普遍的观点^[2，7]。

随着航空航天科技发展，高速流动的数值模拟在飞行器设计验证中起到越来越重要的作用^[8]，CFD领域涌现了多种可压缩流动的处理格式与方法，例如黎曼问题求解器和TVD格式^[9]、用于高精度激波捕捉的WENO格式^[10]，以及计算激波-湍流相互作用的Compact-WENO混合格式^[11]；这些数值方法在理论研究和工程应用领域都取得了可观的成果^[12-13]。不可否认的是，尽管传统CFD方法已经过了多年的发展和完善，在具体问题中要同时准确解析流动小尺度结构和激波间断，实现模拟方法的精度、耗散、稳定性和计算效率等因素的最佳平衡，仍然存在着诸多挑战^[14]。另一方面，工程传热流动问题中，高效换热器普遍具有形状复杂的管腔结构，处理复杂几何边界、并兼顾精度和计算效率也对流动模拟方法提出了更高的要求。

LBM源自于气体动理学理论，对气体运动描述的物理背景清晰，且在处理复杂几何外形、特别是多孔结构时具有独特优势，理论上LBM在模拟可压缩流动和复杂结构的传热流动问题方面拥有巨大潜力。研究者们对传统LBM模型进行了诸多改进，以期突破LBM在计算可压缩流动和热力学流动时的局限性，其基本思想是增加LBM离散模型的自由度和约束条件以还原能量守恒方程。改进方案大致分为两类，一类是额外增加一个能量（内能或总能）分布函数以满足能量方程约束^[15-20]，称之为双分布函数模型。然而，根据气体动理学理论，气体的密度、动量和能量都是气体分子热运动的统计结果，其动力学状态也应当只需由单一的分布函数描述（对于单组分、单原子气体而言），因此引入能量的分布函数必然需要与密度分布函数进行耦合，该耦合关系带有一定的经验性，在增加复杂度的同时也对模型的广泛适用性带来了一定的限制。另一种方案是构造更高阶的单一平衡态分布函数，在平衡态分布函数中引入更多的待定系数，并使用更多的离散速度点，通过增加离散点数量（未知数个数）的方案来满足额外的能量方程约束，典型做法是在每个离散速度方向上使用2个以上的速度大小模态，称之为多层速度模型。与双分布函数模型类似，多层速度模型在能量方程约束条件、分布函数的形式和离散速度点选取上也具有一定的经验性，自Qian^[21]和Alexander等^[22]提出以来，出现了形式多样的方案且各有优缺点。

多层速度模型按平衡态分布函数展开式的构造方法主要分为三类：第一类遵循传统LBM模型的思路，将Maxwell-Boltzmann分布函数泰勒展开至更高阶数^[23-26]。第二类由Shan等^[27]提出，他们发展了Grad^[28-29]的早期工作思路，认为平衡态分布函数可由Hermite多项式级数展开，只需要将展开级数保留足够的截断阶数，即可分别还原到Navier-Stokes方程的质量、动量和能量方程，而密度、动量、能量等宏观统计量均可由Hermite展开式的各阶系数组合表达，并在进一步的工作中阐明传统等温弱可压LBM是Hermite多项式模型的低阶形态，且离散速度是五阶精度Gauss-Hermite积分公式的积分点的事实^[30]。第三类由Xu等^[31]提出，其基本思想是不预设平衡态分布函数的展开形式，而是将离散模型还原质量、动量和能量方程的约束条件视为平衡态分布函数离散点的线性代数方程组，直接求解方程组得出离散分布函数。

由于Maxwell-Boltzmann分布函数的各阶Hermite多项式展开是确定的，Gauss-Hermite积分公式的积分点也具有确定性，使得上述第二类基于Hermite多项式构造的多层速度模型不依赖于经验而具有普适性。然而，由于该模型的数学基础理论较为晦涩，各部分细节的阐述散落于跨度接近20年的多篇文献中，使得科研工作者对该模型进行深入理解存在一定的困难。为了弥补这个缺点，使读者连贯、系统地理解Hermite多项式模型的构造原理和方法，同时也在与其他LBM模型的对比中全面认识各类方案的优劣，我们对Hermite多项式模型重新梳理归纳的同时，也一并列出了其他多层速度模型的构造思想，作为LBM多层速度模型的一个简要总结性工作。

本文的各部分内容为：第1节为Boltzmann方程和LBM的基础理论。由于多层速度模型起源于经典的等温弱可压LBM，且经典LBM模型的很多思想在多层速度模型中仍然适用，故第2节对等温弱可压LBM进行简单介绍。第3节介绍各类多层速度LBM模型（除Hermite多项式模型以外）的基本思想。第4节详细描述Hermite多项式模型的理论和平衡态分布函数的具体形式，以及基于Gauss-Hermite积分公式的离散速度模型及其构造方法。第5节对LBM和传统CFD方法的结合进行了简要介绍，例如LBM有限差分、LBM有限体积和LBM有限元方法。第6节对现有的多层速度模型存在的问题进行总结，并提出未来需要完善和发展的方向。

1 Boltzmann-BGK方程

我们考虑D维（D = 1、2、3）空间中的气体，其中气体分子的位置空间坐标为x，速度为 ${\boldsymbol\xi }$ ，在t时刻的分子密度分布函数为 $ f({{\boldsymbol{x}}},{{\boldsymbol{\xi}} },t) $ 。 $f({{\boldsymbol{x}}},{{\boldsymbol{\xi}} },t){\text{d}}{{\boldsymbol{x}}}{\text{d}}{\boldsymbol\xi }$ 代表t时刻，分布在空间坐标为 ${{\boldsymbol{x}}}～{{\boldsymbol{x}}} + {\text{d}}{{\boldsymbol{x}}}$ 、速度为 ${\boldsymbol\xi }～{\boldsymbol\xi } + {\text{d}}{\boldsymbol\xi }$ 微元内的气体分子质量之和（也就是单个分子质量 ${m_M}$ 与该微元内分子数的乘积）。表征气体微团动力学的宏观统计量，例如密度 $\;\rho $ 、速度 ${{\boldsymbol{u}}}$ 、单位质量的内能e可由 $ f $ 相应表征为：

$ \rho ({{\boldsymbol{x}}},t) = \int {f\;{\text{d}}{\boldsymbol\xi }} $

(1)

$ {{\boldsymbol{u}}}({{\boldsymbol{x}}},t) = \frac{1}{\rho }\int {f{\boldsymbol\xi }\;{\text{d}}{\boldsymbol\xi }} $

(2)

$ e({{\boldsymbol{x}}},t) = \frac{1}{\rho }\int {f{\mkern 1mu} \frac{{|{\boldsymbol\xi } - {{\boldsymbol{u}}}{|^2}}}{2}\;{\text{d}}{\boldsymbol\xi }} $

(3)

以上各式中积分范围为 ${\boldsymbol\xi } \in {( - \infty ,\infty )^D}$ ；由于下文中关于 ${\boldsymbol\xi }$ 的积分范围均与此处相同，故全文中略去积分上下限，且不会导致歧义。定义分子热运动的相对速度 ${{\boldsymbol{v}}} = {\boldsymbol\xi } - {{\boldsymbol{u}}}({{\boldsymbol{x}}},t)$ 。在这些定义的基础上， $ f $ 关于 ${\xi }$ 的二阶矩和三阶矩可表达为：

$ \int {f{\xi _i}{\xi _j}{\text{d}}{\boldsymbol\xi = }\rho {u_i}{u_j} + } \int {f{v_i}{v_j}{\text{d}}{\boldsymbol\xi }} = \rho {u_i}{u_j} - {P_{ij}} $

(4)

$ \begin{split} \frac{1}{2}\int {f{\xi _i}{\xi _j}{\xi _j}{\text{d}}{\boldsymbol\xi}=}&{\rho E{u_i} - } {P_{ij}}{u_j} + \frac{1}{2}\int {f{v_i}{v_j}{v_j}{\text{d}}{\boldsymbol\xi }} \\ =&\rho E{u_i} - {P_{ij}}{u_j} + {q_i} \end{split} $

(5)

其中：

$ {P_{ij}} = - \int {f{v_i}{v_j}{\text{d}}{\boldsymbol\xi }} ,\;\;{q_i} = \frac{1}{2}\int {f{v_i}{v_j}{v_j}{\text{d}}{\boldsymbol\xi }} $

(6)

分别为动量输运项（切应力）和能量扩散项（热流）， $ E = {u^2}/2 + e $ 为单位质量的总能。此处与下文中，对于任意矢量χ有 ${\chi ^2} = {\chi _i}{\chi _i},$ 相同的下标为Einstein约定求和。

根据分子动理学理论和内能e的定义（式（3）），流体的压强p定义为：

$ p = - \frac{1}{D}{P_{ii}} = \frac{1}{D}\int {fv_i^2{\text{d}}{\boldsymbol\xi }} = \frac{2}{D}\rho e $

(7)

能均分定理可表达为：

$ e = \frac{{D{k_B}T}}{{2{m_M}}} = \frac{D}{2}\theta $

(8)

其中，T为气体温度， ${k_B}$ 为Boltzmann常数， $\theta = {k_B}T/{m_M}$ 为无量纲温度。

结合式（7）与式（8）可得到理想气体状态方程：

$ p = \rho \frac{{{k_B}}}{{{m_M}}}T = \rho \theta $

(9)

Boltzmann-BGK方程可表达为：

$ \frac{{\partial f}}{{\partial t}} + {\boldsymbol{\xi} } \cdot {\nabla} f + \frac{{{\boldsymbol{F}}}}{\rho } \cdot {\nabla _{\boldsymbol{\xi}} }f = \varOmega = - \frac{{f - {f^{(0)}}}}{\tau } $

(10)

其中，F为单位体积流体所受的外力， $ \varOmega $ 为碰撞项， $ \tau $ 为粒子因相互碰撞到达平衡态的特征时间（松弛时间），方程的最右端为BGK碰撞模型。 $ {f^{(0)}}({{\boldsymbol{x}}},{\xi },t) $ 为平衡态分布函数，满足Maxwell-Boltzmann分布：

$ {f^{(0)}}({{\boldsymbol{x}}},{\boldsymbol\xi },t) = \frac{\rho }{{{{(2\text{π} \theta )}^{D/2}}}}{{\text{e}}^{ - {{({\boldsymbol\xi } - {u})}^2}/(2\theta )\;}} $

(11)

为了方便之后的分析，我们定义参考尺度 ${L_0}$ 、参考温度 ${T_0}$ 、参考速度 ${V_0} \equiv \sqrt {{k_B}{T_0}/{m_M}} $ 、参考密度 $\;{\rho _0}$ 。对各个量进行无量纲化，则无量纲化的各物理量（以星号表示）为：

$ \begin{split}& V_0^* = 1 \\& {f^*} = fV_0^D/{\rho _0} \\& {{\boldsymbol\xi }^*} = {\boldsymbol\xi }/{V_0} \\& {{{\boldsymbol{x}}}^*} = {{\boldsymbol{x}}}/{L_0} \\& {{{\boldsymbol{u}}}^*} = {{\boldsymbol{u}}}/{V_0} \\& {{{\boldsymbol{v}}}^*} = {{\boldsymbol{v}}}/{V_0} \\& {{{\boldsymbol{F}}}^*} = {{\boldsymbol{F}}}{L_0}/({\rho _0}V_0^2) \\& {t^*} = t{L_0}/{V_0} \\& {\theta ^*} = \theta /V_0^2 \\& {p^*} = p/({\rho _0}V_0^2) \\& {\rho ^*} = \rho /{\rho _0} \end{split} $

(12)

由于无量纲化的Boltzmann-BGK方程和平衡态分布函数分别与式（10）和式（11）形式上完全一致。因此后文无明确标识时，各物理量均为无量纲量，并略去星号标识。对于等温流动，如果取参考温度 ${T_0}$ 为流体本身的温度T，则有无量纲化的 $\theta = 1$ 。由于外力项具有单独的处理方案^[30,32-33]，忽略它不影响对LBM离散模型的分析，下文中将讨论不含外力项的Boltzmann-BGK方程。

对 ${f^{(0)}}$ 分别求 ${\boldsymbol\xi }$ 的零到三阶矩，有：

$ \int {{f^{(0)}}{\text{d}}{\boldsymbol\xi = }\rho } $

(13)

$ \int {{f^{(0)}}{\boldsymbol\xi }{\text{d}}{\boldsymbol\xi }} = \rho {{\boldsymbol{u}}} $

(14)

$ \int {{f^{(0)}}{\xi _i}{\xi _j}{\text{d}}{\boldsymbol\xi }} = \rho {u_i}{u_j} + \rho \theta {\delta _{ij}} $

(15)

$ \begin{split}& \int {{f^{(0)}}{\xi _i}{\xi _j}{\xi _k}{\text{d}}{\boldsymbol\xi }} = \rho {u_i}{u_j}{u_k} +\\&\qquad \quad \rho \theta ({u_i}{\delta _{jk}} + {u_j}{\delta _{ik}} + {u_k}{\delta _{ij}}) \end{split} $

(16)

其中 ${\delta _{ij}}$ 为Kronecker delta函数。在f和 ${f^{(0)}}$ 关于 ${\boldsymbol\xi }$ 各阶矩的基础上，我们对Boltzmann-BGK方程的两端分别乘以1、 ${\xi _i}$ 、 ${\xi _j}{\xi _j}$ 并积分，可得到质量、动量和能量守恒方程：

$ \frac{{\partial \rho }}{{\partial t}} + \frac{{\partial (\rho {u_i})}}{{\partial {x_i}}} = 0 $

(17)

$ \frac{{\partial \rho {u_i}}}{{\partial t}} + \frac{{\partial (\rho {u_i}{u_j})}}{{\partial {x_j}}} = \frac{{\partial {P_{ij}}}}{{\partial {x_j}}} $

(18)

$ \frac{{\partial (\rho E)}}{{\partial t}} + \frac{{\partial (\rho {u_i}E)}}{{\partial {x_i}}} = \frac{{\partial ({P_{ij}}{u_j})}}{{\partial {x_i}}} - \frac{{\partial {q_i}}}{{\partial {x_i}}} $

(19)

可以看到，式（17）～式（19）已经十分接近Navier-Stokes方程。LBM中选定f的离散模型，将决定 ${P_{ij}}$ 、 ${q_i}$ 的表达和该模型可适用的范围。

2 等温弱可压LBM模型

等温弱可压LBM模型是最早提出的LBM模型，并在不可压流动计算中得到广泛应用。严格意义上，包括气体、液体在内的所有流体都具有可压缩性，例如液体的压缩性可用Tait方程描述^[34]。传统求解压力泊松方程的不可压流计算方法将导致无限大声速，也就是压力扰动的无限速度传播，这与物理事实是相悖的。采用弱可压方法计算“不可压流动”理论上更加接近物理本质。本节对等温弱可压LBM模型进行简要陈述，这也将成为多层速度LBM模型修正和拓展的基础。

2.1 速度空间离散模型

LBM的核心是将连续的Boltzmann方程离散化，即用速度空间离散的分布函数 ${f_a}({{\boldsymbol{x}}},t)$ 和 $f_a^{(0)}({{\boldsymbol{x}}},t)$ 分别代替连续的分布函数 $f({{\boldsymbol{x}}},{\boldsymbol\xi },t)$ 和 ${f^{(0)}}({{\boldsymbol{x}}},{\boldsymbol\xi },t),$ 其中a = 1,2,3,···,d为离散点序号。为此形成了两种思路：泰勒展开法（也称之为马赫展开法）和Hermite多项式展开法。此处介绍泰勒展开法，而Hermite多项式展开法将在第4节Hermite多项式模型部分详细介绍。

将平衡态分布函数 ${f}^{(0)}({{\boldsymbol{x}}},\boldsymbol\xi,t)$ 泰勒展开到u的二阶项：

$ \begin{split} {f}^{(0)} =& \frac{\rho }{{(2\text{π} \theta )}^{D/2}}{\text{e}}^{-\text{}\text{}\text{}\text{}\frac{{(\xi -u)}^{2}}{2\theta }}\\ =& \frac{\rho }{{(2\text{π} \theta )}^{D/2}}{\text{e}}^{-\frac{{\xi}^{2}}{2\theta }}\mathrm{exp}\left(\frac{\xi\cdot {\boldsymbol{u}}}{\theta }-\frac{{u}^{2}}{2\theta }\right)\\ =& \frac{\rho }{{(2\text{π} \theta )}^{D/2}}{\text{e}}^{-\frac{{\xi}^{2}}{2\theta }}\left[1+\frac{\xi\cdot {\boldsymbol{u}}}{\theta }+\frac{{\left(\xi\cdot {\boldsymbol{u}}\right)}^{2}}{2{\theta }^{2}}-\frac{{{u}}^{2}}{2\theta }\right] \end{split}$

(20)

将 ${f}^{(0)}({{\boldsymbol{x}}},\;\boldsymbol\xi,\;t)$ 在速度空间的系列离散点 ${{\boldsymbol\xi }_a}$ 代入式（20），得到离散的平衡态分布函数：

$ f_a^{(0)} = \rho {w_a}\left[ {1 + \frac{{{{\boldsymbol\xi }_a} \cdot {{\boldsymbol{u}}}}}{\theta } + \frac{{{{\left( {{{\boldsymbol\xi }_a} \cdot {{\boldsymbol{u}}}} \right)}^2}}}{{2{\theta ^2}}} - \frac{{{{{\boldsymbol{u}}}^2}}}{{2\theta }}} \right] $

(21)

其中， ${w_a} = {(2\text{π} \theta )^{ - D/2}}{{{\rm{exp}}}}\left[ - \boldsymbol\xi _a^2/(2\theta )\right]$ 为对应于离散速度 $ {{{\boldsymbol{\xi }}}_a} $ 的权重。 $ {{{\boldsymbol{\xi}} }_a} $ 和 $ {w_a} $ 需要满足的条件为：离散平衡态分布函数 $f_a^{(0)}$ 对 $ {{\boldsymbol{\xi}} } $ 的零到二阶矩（分别为质量、动量和动量输运量）与连续平衡态分布函数 $ {f^{(0)}} $ 对 $ {{\boldsymbol{\xi}} } $ 的零到二阶矩相等，即：

$ \rho = \sum\limits_{a = 1}^d {f_a^{(0)} = \int {{f^{(0)}}{\mkern 1mu} {\text{d}}{{\boldsymbol{\xi}} }} } $

(22)

$ \rho {{\boldsymbol{u}}} = \sum\limits_{a = 1}^d {f_a^{(0)}{{{\boldsymbol{\xi}} }_a} = \int {{f^{(0)}}{{\boldsymbol{\xi}} }{\mkern 1mu} {\text{d}}{{\boldsymbol{\xi}} }} } $

(23)

$ \rho {u_i}{u_j} + \rho \theta {\delta _{ij}} = \sum\limits_{a = 1}^d {f_a^{(0)}{\xi _{ai}}{\xi _{aj}} = \int {{f^{(0)}}{\xi _i}{\xi _j}{\mkern 1mu} {\text{d}}{{\boldsymbol{\xi}} }} } $

(24)

将式（21）代入式（22）～式（24），可得到关于 $ {{{\boldsymbol{\xi}} }_a} $ 和 $ {w_a} $ 的约束条件^[35-37]：

$ \sum\limits_{a = 1}^d {{w_a} = 1} $

(25)

$ \sum\limits_{a = 1}^d {{w_a}{\xi _{ai}}{\xi _{aj}} = \theta {\delta _{ij}}} $

(26)

$ \sum\limits_{a = 1}^d {{w_a}{\xi _{ai}}{\xi _{aj}}{\xi _{ak}}{\xi _{al}} = {\theta ^2}({\delta _{ij}}} {\delta _{kl}} + {\delta _{ik}}{\delta _{jl}} + {\delta _{il}}{\delta _{jk}}) $

(27)

此时我们根据经验预设对称性的离散速度模型，例如图1所示的D1Q3（ $d = 3$ 、 ${\xi _1} = 0$ 、 ${\xi _2} = - {\xi _3}= c{V_0}$ 、 ${w_2} = {w_3}$ ）和D2Q9模型等，并将离散速度模型和 $\theta = 1$ 代入式（25）～式（27）求解。对于D1Q3可解得 $ c = \sqrt{3} $ 、 $ {w}_{1} = 2/3 $、 ${w_2} = 1/6$ 。类似的方法应用到二维和三维情况预设的离散速度模型中，同样可解出相应的 $c$ 和权重系数，进而得到具体的D2Q9、D3Q19、D3Q27等模型。

2.2 时间和空间离散

在获得离散速度 ${{{\boldsymbol{\xi}} }_a}$ 和权重 $ {w_a} $ 后，将 ${{\boldsymbol{\xi}}}_{a}、$ $ {w_a} $ 和 $ f_a^{(0)} $ 代入Boltzmann-BGK方程（式（10）），得到关于 ${f_a} = f({{\boldsymbol{x}}},{{\boldsymbol\xi }_a},t)$ 的离散方程；并将方程的左端视为关于 $ {f_a} $ 的物质导数，沿各离散速度 ${{\boldsymbol\xi }_a}$ 所代表的特征线方向对f进行时间积分，有：

$ \begin{split}& {f_a}({{\boldsymbol{x}}} + {{\boldsymbol\xi }_a}\Delta t,t + \Delta t) - {f_a}({{\boldsymbol{x}}},t) = \int_t^{t + \Delta t} {{\varOmega _a}({{\boldsymbol{x}}},t){\text{d}}t} \end{split} $

(28)

使用梯形积分公式达到二阶时间精度，式（28）可表达为：

$ \begin{split}&\quad {f_a}({{\boldsymbol{x}}} + {{\boldsymbol\xi }_a}\Delta t,t + \Delta t) - {f_a}({{\boldsymbol{x}}},t) \\& = \frac{{\Delta t}}{2}[{\varOmega _a}({{\boldsymbol{x}}},t) + {\varOmega _a}({{\boldsymbol{x}}} + {{\boldsymbol\xi }_a}\Delta t,t + \Delta t)] \\& = - \frac{{\Delta t}}{{2\tau }}[{f_a}({{\boldsymbol{x}}},t) - f_a^{(0)}({{\boldsymbol{x}}},t) + \\&\quad {f_a}({{\boldsymbol{x}}} + {{\boldsymbol\xi }_a}\Delta t,t + \Delta t) - f_a^{(0)}({{\boldsymbol{x}}} + {{\boldsymbol\xi }_a}\Delta t,t + \Delta t)] \end{split} $

(29)

为使得式（29）变为显式方法，定义中间时刻的分布函数^[16]：

$ {\overline f _a} = {f_a} + \frac{{\Delta t}}{{2\tau }}({f_a} - f_a^{(0)}) $

(30)

则式（29）变为：

$ \begin{split}& {\overline f _a}({{\boldsymbol{x}}} + {{\boldsymbol\xi }_a}\Delta t,t + \Delta t) - {\overline f _a}({{\boldsymbol{x}}},t) = - \frac{{\Delta t}}{{\bar \tau }}[{\overline f _a}({{\boldsymbol{x}}},t) - f_a^{(0)}({{\boldsymbol{x}}},t)] \end{split} $

(31)

其中：

$ \bar \tau = \tau + \Delta t/2 $

(32)

使用Chapman-Enskog展开分析，并对离散速度模型和平衡态分布函数添加额外约束条件：

$ \begin{split} \sum\limits_{a = 1}^d {f_a^{(0)}{\xi _{ai}}{\xi _{aj}}{\xi _{ak}} =}& {\int {{f^{(0)}}{\xi _i}{\xi _j}{\mkern 1mu} {\xi _k}{\mkern 1mu} {\text{d}}{\boldsymbol\xi }} - \rho {u_i}{u_j}{u_k}} \\= & \rho \theta ({u_i}{\delta _{jk}} + {u_j}{\delta _{ik}} + {u_k}{\delta _{ij}}) \end{split} $

(33)

可证明式（31）恢复了Navier-Stokes方程的质量和动量方程^[33,36]，且流体的运动学黏性系数为：

$ \nu = \frac{p}{\rho }\Bigg(\bar \tau - \frac{{\Delta t}}{2}\Bigg) $

(34)

需要注意的是，由于约束条件式（33）成立的前提是u的三阶量可忽略，因此Chapman-Enskog展开也只能在忽略 ${O}({{{\boldsymbol{u}}}^3})$ 量的情况下才能获得准确的黏性应力项，因此本节介绍的方法只适用于 ${{\boldsymbol{u}}} \to 0$ 的情形，即等温低速弱可压流动。

常用的D1Q3和D2Q9等离散速度模型中， ${{\boldsymbol\xi }_a}$ 的各分量与单位速度 ${V_0}$ 之间相差一个倍数c，此时将所有的速度 ${{\boldsymbol\xi }_a}$ 和 $ {{\boldsymbol{u}}} $ 都乘以倍数 $ {c_s} = 1/c $ 得到归一化速度： ${{\boldsymbol\xi }_a}' = {\boldsymbol{c}_a} \equiv {c_s}{{\boldsymbol\xi }_a}$ ， $ {{\boldsymbol{u}}'} \equiv {c_s}{{\boldsymbol{u}}} $ ，此时 $ {{{\boldsymbol{c}}}_a} $ 各分量只有0或±1三种形态，在编程实现中更加便利。式（21）用归一化速度表达为：

$ f_a^{(0)} = \rho {w_a}\left[ {1 + \frac{{{\boldsymbol{c}_a} \cdot {{\boldsymbol{u}}'}}}{{c_s^2\theta }} + \frac{{{{\left( {{\boldsymbol{c}_a} \cdot {{\boldsymbol{u}}'}} \right)}^2}}}{{2c_s^4{\theta ^2}}} - \frac{{{{\boldsymbol{u}}}{{'}^2}}}{{2c_s^2\theta }}} \right] $

(35)

相应于归一化速度，我们定义归一化时间 $t' = t/{c_s}$ ，在LBM计算中，一般取无量纲化的网格宽度 $\Delta x = 1$ ；此时对归一化的时间步 $ \Delta t' = 1 $ ，实际为 $ \Delta t = {c_s} $ ，于是 $ t $ 时刻坐标为x的网格上速度为 ${{\boldsymbol\xi }_a}$ 的粒子在经过 $ \Delta t $ 时间的传播后，恰好也落在相邻网格 ${{\boldsymbol{x}}} {+} {\Delta} {{\boldsymbol{x}}}$ 上。另外，由于压强为分布函数的速度二阶矩，归一化的压强为 $p' = c_s^2p$ 。由于归一化的各物理量表达和使用上都更为简洁，在各类文献最终的LBM模型都略去表示归一化的符号，将平衡态分布函数离散式表示为：

$ {f}_{a}^{（0）} = \rho {w}_{a}\left[1+\frac{\boldsymbol{c}_{a}\cdot {\boldsymbol{u}}}{{c}_{s}^{2}}+\frac{{\left(\boldsymbol{c}_{a}\cdot {\boldsymbol{u}}\right)}^{2}}{2{c}_{s}^{4}}-\frac{{{\boldsymbol{u}}}^{2}}{2{c}_{s}^{2}}\right] $

(36)

相应的无量纲时间和空间离散式也写为：

$ \begin{split}& {\bar f _a}({{\boldsymbol{x}}} + {\boldsymbol{c}_a},t + 1) - {\bar f _a}({{\boldsymbol{x}}},t) = - \frac{1}{{\bar \tau }}\big[{\bar f _a}({{\boldsymbol{x}}},t) - f_a^{(0)}({{\boldsymbol{x}}},t)\big] \end{split} $

(37)

其中：

$ \bar \tau = \tau + 1/2 $

(38)

相应的无量纲黏性系数为：

$ \upsilon = c_s^2\Bigg(\bar \tau - \frac{1}{2}\Bigg) $

(39)

${f}_{a}^{(0)}$ 中的密度 $\rho $ 和速度 ${{\boldsymbol{u}}}$ 可按式（1）～式（2）的离散形式计算：

$ \rho = \sum\limits_{a = 1}^d {{f_a}} $

(40)

$ {{\boldsymbol{u}} = }\frac{1}{\rho }\sum\limits_{a = 1}^d {{f_a}{{{\boldsymbol{c}}}_a}} $

(41)

式（40）～式（41）的合理性将在第4节详细说明。

3 多层速度LBM模型

从第2节可以看到，等温弱可压流动LBM模型，其离散速度在每个方向都只有一个值（分布函数信息只传播到相邻的网格点），因此称之为单层速度模型。推导该模型时，只约束了离散模型的质量、动量和动量输运量，故无法准确还原能量方程。受此启发，业内学者开始对离散模型施加能量相关（ $ f_{}^{(0)} $ 的 ${\boldsymbol\xi }$ 三阶以上矩）的约束条件，伴随而来的是更高阶的 $ f_{}^{(0)} $ 展开式和更多的离散速度点，其中普遍做法是在各离散速度方向上采用2个以上的速度值（模态），因此这些模型称之为多层速度模型。本节对几类较为典型的多层速度模型进行简要陈述。

3.1 早期多层速度模型

Qian^[21]和Alexander等^[22]分别在1993年提出了最早的多层速度模型。他们均指出， ${f^{(0)}}$ 的展开形式需要添加u的三阶量，其中Alexander等给出的 ${f^{(0)}}$ 的离散形式为：

$ \begin{split} f_{pa}^{(0)} =& \rho \big[{A_p} + {B_p}{{\boldsymbol\xi }_a} \cdot {{\boldsymbol{u}} + }{C_p}{({{\boldsymbol\xi }_a} \cdot {{\boldsymbol{u}}})^2} + {D_p}{u^2} { + } \\& {E_p}{({{\boldsymbol\xi }_a} \cdot {{\boldsymbol{u}}})^3} + {F_p}{u^2}{{\boldsymbol\xi }_a} \cdot {{\boldsymbol{u}}}\big] \end{split} $

(42)

其中，p = 0,1,2为预设离散速度模型的网格层数，p = 0为中心网格，p = 1，2为中心周围第1、第2层网格；a = 1, 2, 3,···m为该层网格上的离散点编号。相应的二维离散速度模型D2Q13如图2～图3所示。

对 ${f^{(0)}}$ 离散式及离散速度模型的约束条件除了式（22）～式（24）以外（对所有的离散点求和），还需增加用离散速度表征的热流项：

$ {q_i} = \frac{1}{2}\sum\limits_{p,a}^{} {{f_{pa}}} {({\xi _{paj}} - {u_j})^2}({\xi _{pai}} - {u_i}) $

(43)

在Chapam-Enskog展开中， ${q_i}$ 的零阶近似为0，一阶近似与温度梯度成正比。在这些约束条件下可解出 ${A_p}～{F_p}$ 的表达式（具体的求解方法可参考Wolf-Gladrow^[35]附录6.4）。

Chen等^[23]认为，Qian和Alexander等提出的模型，在还原Navier-Stokes的能量方程后会额外多出非线性误差项 ${\partial _i}{\partial _j}(\rho e{u_i}{u_j})$ 、 ${\partial _i}{\partial _j}(\rho {u_i}{u_j})$ 等。他们按照Frish^[1]的思路，定义离散速度张量：

$ \boldsymbol{T}_{pk{\alpha _1}{\alpha _2}{\alpha _2} \cdots {\alpha _n}}^{(n)} = \sum\limits_{i = 1}^m {{\xi _{pki{\alpha _1}}}} {\xi _{pki{\alpha _2}}}{\xi _{pki{\alpha _3}}} \cdots {\xi _{pki{\alpha _n}}} $

(44)

其中：（p，k）为组编号，代表该组离散速度的大小均为 $k\sqrt p $ ，且共有m个对称的离散速度点； ${\alpha _1}～{\alpha _n}$ 为n阶张量的坐标下标；该速度张量需要保证二、四、六阶均为各向同性的（由于离散速度的坐标对称性奇数阶张量自动为0）。 ${f^{(0)}}$ 的离散式需添加速度的四阶项：

$ \begin{split} f_{pka}^{(0)} =& {A_{pk}} + {M_{pk}}{{\boldsymbol\xi }_{pka}} \cdot {{\boldsymbol{u}} + }{J_{pk}}{({{\boldsymbol\xi }_{pka}} \cdot {{\boldsymbol{u}}})^2} + \\& {G_{pk}}{u^2}{ + }{H_{pk}}{({{\boldsymbol\xi }_{pka}} \cdot {{\boldsymbol{u}}})^3} + {Q_{pk}}{u^2}{{\boldsymbol\xi }_{pka}} \cdot {{\boldsymbol{u}}} { + } \\& {R_{pk}}{u^2}{({{\boldsymbol\xi }_{pka}} \cdot {{\boldsymbol{u}}})^2} + {S _{pk}}{u^4} \end{split} $

(45)

对 ${f^{(0)}}$ 的离散式和离散速度约束条件，除了 $f_{pka}^{(0)}$ 的速度零到二阶矩满足式（22）～式（24）以外，其速度三阶矩需满足：

$ \begin{split} \sum\limits_{p,k,a}^{} {f_{pka}^{(0)}{\xi _{ai}}{\xi _{aj}}{\xi _{ak}}} =& \int {{f^{(0)}}{\xi _i}{\xi _j}{\mkern 1mu} {\xi _k}{\mkern 1mu} {\text{d}}{\boldsymbol\xi }} \\ =& \rho {u_i}{u_j}{u_k} + \rho \theta ({u_i}{\delta _{jk}} + {u_j}{\delta _{ik}} + {u_k}{\delta _{ij}}) \end{split} $

(46)

此时离散模型与连续的 $ {f^{(0)}} $ 速度三阶矩一致，这对于还原能量输运项至关重要。同时， $f_{pka}^{(0)}$ 的速度四阶矩需满足：

$ \begin{split}& \sum\limits_{p,k,a}^{} {f_{pka}^{(0)}\xi _a^2{\xi _{ai}}{\xi _{aj}}} = \rho \left[ {\frac{{4(D + 2)}}{{{D^2}}}} \right.{e^2}{\delta _{ij}} +\\&\qquad \frac{2}{D}e{u^2}{\delta _{ij}} + \left. {\frac{{2(D + 4)}}{D}e{u_i}{u_j} + {u^2}{u_i}{u_j}} \right] \end{split} $

(47)

在以上约束条件下，预设一维多层速度模型1D5V，其离散点为 ${\xi }_{a} = \{0，\pm c，\pm 2c\}$ 。二维离散速度为一、二层网格除中心点外的全部速度点，如图4所示，共16个，称之为2D16V；类似的，三维多层速度模型为3D40V。

3.2 Watari-Tsutahara模型

Watari和Tsutahara的计算试验发现^[24-25]，Chen等的模型在计算热流和可压缩问题时仍然不能令人满意。例如在计算平板Couette流动时，一旦松弛时间 $\tau < 0.2$ ，无论上下平板的温度是否一致，计算结果与理论解之间都存在较明显的差异，这说明Chen等的模型只在较窄的黏性系数范围内成立。他们沿用了Chen等提出的约束条件，对 ${f^{(0)}}$ 的离散形式作了进一步改进，提出了两种方案：一种是基于泰勒展开的思路，将Maxwell-Boltzmann分布展开到u的四阶项：

$ \begin{split} f_{pka}^{(0)} =& \rho {F_{pk}}\left[ {\left( {1 - \frac{D}{{4e}}{u^2} + \frac{{{D^2}}}{{32{e^2}}}{u^4}} \right)} \right. \\& + \frac{D}{{2e}}\left( {1 - \frac{D}{{4e}}{u^2}} \right){{\boldsymbol\xi }_{pka}} \cdot {{\boldsymbol{u}} + }\frac{{{D^2}}}{{8{e^2}}}\left( {1 - \frac{D}{{4e}}{u^2}} \right){({{\boldsymbol\xi }_{pka}} \cdot {{\boldsymbol{u}}})^2} \\& + \frac{{{D^3}}}{{48{e^3}}}{({{\boldsymbol\xi }_{pka}} \cdot {{\boldsymbol{u}}})^3} + \left. {\frac{{{D^4}}}{{384{e^4}}}{{({{\boldsymbol\xi }_{pka}} \cdot {{\boldsymbol{u}}})}^4}} \right] \end{split} $

(48)

另一种是在Chen等的模型，即式（45）中再添加 ${\boldsymbol\xi }$ 的4次项：

$ \begin{split} f_{pka}^{(0)} =& \rho [{A_{pk}} + {M_{pk}}{{\boldsymbol\xi }_{pka}} \cdot {{\boldsymbol{u}} + }{J_{pk}}{({{\boldsymbol\xi }_{pka}} \cdot {{\boldsymbol{u}}})^2} \\& + {G_{pk}}{u^2}{ + }{H_{pk}}{({{\boldsymbol\xi }_{pka}} \cdot {{\boldsymbol{u}}})^3} + {Q_{pk}}{u^2}{{\boldsymbol\xi }_{pka}} \cdot {{\boldsymbol{u}}} \\& + {R_{pk}}{u^2}{({{\boldsymbol\xi }_{pka}} \cdot {{\boldsymbol{u}}})^2} + {S _{pk}}{u^4} + {T_{pk}}{({{\boldsymbol\xi }_{pka}} \cdot {{\boldsymbol{u}}})^4}] \end{split} $

(49)

需要说明的是，在Watari和Tsutahara^[25]的原文中，将 $|{{\boldsymbol\xi }_{pka}}|$ 相同的离散速度点组编号简化为单序号 $k$ ，但其离散速度模型中依然以2位数表示 $k$ 。本文中为了与Chen^[23]的表达方法对照，依然沿用了（p，k）来表示组编号。这两种模型采用的离散速度模型具有4层网格，其二维模型如图5所示。

3.3 比热可变模型

本小节之前介绍的几类多层速度模型，事实上都是基于单原子气体动理学理论，这些模型会导致还原的气体比热比为固定值，不适用于双原子和多原子分子气体。为了弥补这一缺陷，Kataoka和Tsutahara提出^[26]，以 ${E_b} = (b{k_B}T/{m_M} + {u^2})/2$ 替换单原子分子的总能E，其中常数b为气体分子的平动和转动、振动的总自由度，与比热比 $\gamma $ 的关系为 $\gamma = b/(b + 2)$ ，并引入 ${\eta ^2}$ 用于描述平动以外自由度的内能，此时温度T用离散速度点（ ${{\boldsymbol\xi }_a} $ ， ${\eta _a} $ ）定义为：

$ \rho (b{k_B}T/{m_M} + {u^2}) = \sum\limits_{a = 1}^d {{f_a}(\xi _a^2 + \eta _a^2)} $

(50)

对离散速度和 $f_a^{(0)}$ 的约束条件，除了式（22）～式（24），其余的为：

$ \rho (b{k_B}T/{m_M} + {u^2}) = \sum\limits_{a = 1}^d {f_a^{(0)}(\xi _a^2 + \eta _a^2)} $

(51)

$ \rho {u_i}{u_j}{u_k} + \rho \theta ({u_i}{\delta _{jk}} + {u_j}{\delta _{ik}} + {u_k}{\delta _{ij}}) = \sum\limits_{a = 1}^d {f_a^{(0)}{\xi _{ai}}{\xi _{aj}}{\xi _{ak}}} $

(52)

$ \rho [(b + 2){k_B}T/{m_M} + {u^2}]{u_i} = \sum\limits_{a = 1}^d {f_a^{(0)}(\xi _a^2 + \eta _a^2)} {\xi _{ai}} $

(53)

$ \begin{split}& \rho \Big\{ (b + 2){({k_B}T/{m_M})^2}{\delta _{ij}} + \big[(b + 4){u_i}{u_j} + {u^2}{\delta _{ij}}\big]{k_B}T/{m_M} \\&\qquad + {u^2}{u_i}{u_j}\Big\} = \sum\limits_{a = 1}^d {f_a^{(0)}(\xi _a^2 + \eta _a^2)} {\xi _{ai}}{\xi _{aj}} \\[-12pt] \end{split}$

(54)

Kataoka和Tsutahara给出了二维情况下一个16点离散速度模型，如图6所示。而 ${\eta _a}$ 为：

$ {\eta }_{a}/{c}_{s} = \Bigg\{\begin{array}{l}5/2,\;\;\;1\leqslant a\leqslant 4\\ 0,\;\;\;\;\;\;\;5\leqslant a\leqslant 16\end{array} $

(55)

其 $ f_{}^{(0)} $ 离散展开式为：

$ \begin{split} f_a^{(0)} =& \rho \big[{A_{0a}} + {A_{1a}}T + {A_{2a}}{T^2} + ({A_{3a}} + {A_{4a}}T)u_a^2 + \\& {A_{5a}}{(u_a^2)^2} + ({B_{0a}} + {B_{1a}}T + {B_{2a}}{u^2}){{\boldsymbol\xi }_{a}} \cdot {{\boldsymbol{u}}} + \\& ({D_{0a}} + {D_{1a}}T + {D_{2a}}{u^2}){({{\boldsymbol\xi }_{a}} \cdot {{\boldsymbol{u}}})^2} + {E_a}{({{\boldsymbol\xi }_{a}} \cdot {{\boldsymbol{u}}})^3}\big] \end{split} $

(56)

根据约束条件可求出 $ f_a^{(0)} $ 中的所有待定系数。

Xu等^[31,38-39]则在Kataoka-Tsutahara模型的基础上，发展了离散Boltzmann方法（DBM）。传统多层速度模型往往先预设含待定系数的平衡态分布函数展开形式，根据约束条件求得待定系数，再将离散速度点代入 $ f_{}^{(0)} $ 展开式得到 $ f_a^{(0)} $ ；DBM不再预设 $ f_a^{(0)} $ 的形式，而是直接根据约束条件和离散速度模型解出平衡态分布函数的离散点 $ f_a^{(0)} $ 。其具体思路是：将所有7个约束条件式（22）～式（24）、式（51）～式（54）中的向量方程按坐标分量分别列出（例如D = 2时式（53）将被列为2个方程），并统计所有的方程个数I；此时将这I个约束方程，视为以I个平衡态分布函数的离散点 $ f_a^{(0)} $ 为未知数的方程，因此离散速度模型中 ${{\boldsymbol\xi }_a}$ 的数量为I (d = I)。定义未知数向量：

$ {{{\boldsymbol{f}}}^{(0)}} = {(f_1^{(0)},f_2^{(0)},f_3^{(0)}, \cdots ,f_I^{(0)})^{\text{T}}} $

(57)

将约束条件的等式左边部分视为方程组的右端项向量 ${\boldsymbol{M}} = {(\rho ,\rho {u}_{x},\rho {u}_{y},\cdots )}^{\text{T}}，$ 则约束条件可写为线性方程组：

$ C{{\boldsymbol{f}}}^{(0)} = {{\boldsymbol{M}}} $

(58)

求解方程组即可得到分布函数的离散形式 ${{{\boldsymbol{f}}}^{(0)}} = {C^{ - 1}}{{\boldsymbol{M}}}$ 。对于二维情况，I = 16，为此对 ${{\boldsymbol\xi }_a}$ 而言可使用多种16点离散速度模型，如图7所示。对于平动以外的自由度 ${\eta _a}$ ，其选取的基本原则是系数矩阵 $C$ 可逆，以保证方程（58）有解；且根据能均分定理，需要 ${\eta _a}$ 与 $\sqrt {(b - D)\theta } $ 在同一数量级^[40]。该方法已在可压缩流动、燃烧、流体界面不稳定性等多种非平衡态、多相流复杂系统模拟中得到应用，并取得较满意的效果^[41-42]。

4 Hermite多项式模型 4.1 速度空间展开

从第3节的介绍可以看出，大部分多层速度模型还是沿袭了单层速度模型“预设平衡态分布函数离散形式—预设离散速度模型—待定系数法求解”的思路。这种思路可以针对每一类特定问题都构造相对优化的多层速度模型，但经验性较强，对于开始接触多层速度模型的学者而言，面对种类繁多的离散模型、平衡态分布函数形式和约束条件，如何选择最优方案往往无所适从；即便获得了对某一问题的优化组合，当流动参数改变后，计算结果不一定同样令人满意，无法保证模型的普适性。

Shan等^[27,30]在Grad^[28-29]工作的基础上，发展了Hermite多项式模型，该模型基于Hermite多项式级数展开的数学性质，使平衡态分布函数的展开形式和离散速度模型都不再依赖于经验而具有确定性，且可以证明单层速度模型事实上是多层速度模型的低阶形态。这正是本节将详述的内容。

参照Grad的定义，D维空间的n阶Hermite多项式定义为n阶张量：

$ \boldsymbol{H}_{{i_1}{i_2} \cdots {i_n}}^{(n)}({\boldsymbol\xi }) = \frac{{{{( - 1)}^n}}}{{\omega ({\boldsymbol\xi })}}\nabla _{{i_1}{i_2} \cdots {i_n}}^n\omega ({\boldsymbol\xi }) $

(59)

其中：

$ \nabla _{{i_1}{i_2} \cdots {i_n}}^n = \frac{\partial }{{\partial {\xi _{{i_1}}}}}\frac{\partial }{{\partial {\xi _{{i_2}}}}} \cdots \frac{\partial }{{\partial {\xi _{{i_n}}}}} $

(60)

$ \omega ({\boldsymbol\xi }) = \frac{1}{{{{(2\text{π} )}^{D/2}}}}\exp \left( { - \frac{{{\boldsymbol\xi } \cdot {\boldsymbol\xi }}}{2}} \right) $

(61)

不同阶数之间的Hermite多项式之间有递推关系：

$ {\xi _i}\boldsymbol{H}_{{i_1}{i_2} \cdots {i_n}}^{(n)} = \boldsymbol{H}_{{i_1}{i_2} \cdots {i_n}i}^{(n + 1)} + \sum\limits_{k = 1}^n {{\delta _{{i_k}i}}\boldsymbol{H}_{{i_1}{i_2} \cdots {i_{k - 1}}{i_{k + 1}} \cdots {i_n}}^{(n - 1)}} $

(62)

以及正交关系：

$ \int {\omega ({\boldsymbol\xi })\boldsymbol{H}_{\alpha }^{(m)}({\boldsymbol\xi })\boldsymbol{H}_{\beta }^{(n)}({\boldsymbol\xi }){\text{d}}\boldsymbol\xi } = \prod\limits_{i = 1}^D {{n_i}!} \;{\delta _{mn}}{\delta _{{\alpha \beta }}} $

(63)

其中， ${\delta _{{\alpha \beta }}}$ 只有在序列 $\;{{\boldsymbol{\beta}} } = ({\beta _1},{\beta _2},{\beta _3}, \cdots {\beta _n})$ 为序列 ${\text{α}} = ({\alpha _1},{\alpha _2},{\alpha _3}, \cdots {\alpha _m})$ 的任意全排列时为1，否则为0； ${n_i}$ 为第i个下标 ${\alpha _i}$ 在整个排列 ${\text{α}}$ 中重复的次数。

根据Hermite多项式的定义和递推关系，零到四阶Hermite多项式为：

$ {H^{(0)}}({{\boldsymbol{\xi}} }) = 1 $

(64)

$ H_i^{(1)}({{\boldsymbol{\xi}} }) = {\xi _i} $

(65)

$ H_{ij}^{(2)}({{\boldsymbol{\xi}} }) = {\xi _i}{\xi _j} - {\delta _{ij}} $

(66)

$ H_{ijk}^{(3)}({{\boldsymbol{\xi}} }) = {\xi _i}{\xi _j}{\xi _k} - {\xi _i}{\delta _{jk}} - {\xi _j}{\delta _{ik}} - {\xi _k}{\delta _{ij}} $

(67)

$ \begin{split} H_{ijkl}^{(4)}({{\boldsymbol{\xi}} }) =& {\xi _i}{\xi _j}{\xi _k}{\xi _l} - ({\xi _i}{\xi _j}{\delta _{kl}} + {\xi _i}{\xi _k}{\delta _{jl}} + \\& {\xi _i}{\xi _l}{\delta _{jk}} + {\xi _j}{\xi _k}{\delta _{il}} + {\xi _j}{\xi _l}{\delta _{ik}} + {\xi _k}{\xi _l}{\delta _{ij}}) + \\& {\delta _{ij}}{\delta _{kl}} + {\delta _{ik}}{\delta _{jl}} + {\delta _{il}}{\delta _{jk}} \end{split} $

(68)

由于Hermite多项式序列为一组完备基，对于任意在 $ {(-\infty ，+\infty )}^{D} $ 区间对权函数 $\omega ({{\boldsymbol{\xi }}})$ 平方可积的函数 $f({{\boldsymbol{\xi }}})$ ，可以用Hermite多项式级数展开为：

$ \begin{split} f({{\boldsymbol{\xi}} }) =& \omega ({{\boldsymbol{\xi}} })\sum\limits_{n = 0}^\infty {\frac{1}{{n!}}{{\boldsymbol{a}}}_{}^{(n)} \cdot {{\boldsymbol{H}}}_{}^{(n)}({{{\boldsymbol{\xi}}} })} = \omega ({{\boldsymbol{\xi}} })\sum\limits_{n = 0}^\infty {\frac{1}{{n!}}\sum\limits_{i} {a_{i}^{(n)}H_{i}^{(n)}({{\boldsymbol{\xi}} })} } \end{split} $

(69)

$ a_{i}^{(n)} = \int {f({{\boldsymbol{\xi}} })} H_{i}^{(n)}({{\boldsymbol{\xi}} }){\text{d}}{{\boldsymbol{\xi}} } $

(70)

在一维情况下，n阶Hermite多项式具有唯一表达式；多维n阶Hermite多项式则是由下标排列 ${{\boldsymbol{i}}} = ({i_1},{i_2},{i_{3,}} \cdots {i_n})$ 确定的D维多项式构成，因此 $ {{\boldsymbol{a}}}_{}^{(n)} = \{ a_{i}^{(n)}\} $ 和 $ {{\boldsymbol{H}}}_{}^{(n)}({\xi }) = \{ H_{i}^{(n)}({\xi })\} $ 为 ${{\boldsymbol{i}}}$ 的全排列决定的对称张量集合。

容易证明， $f({{\boldsymbol{\xi}} })$ 关于速度 ${\boldsymbol\xi }$ 任意N阶矩的积分因子 ${p_{_N}}({{\boldsymbol{\xi}} }) = {\xi _{{i_1}}}{\xi _{{i_2}}}{\xi _{{i_3}}} \cdots {\xi _{{i_N}}}$ 可由不高于N阶的Hermite多项式线性表达出^[43]，即：

$ {p_{_N}}({{\boldsymbol{\xi}} }) = \sum\limits_{n = 0}^N {{{\boldsymbol{b}}}_{}^{(n)} \cdot {{\boldsymbol{H}}}_{}^{(n)}({{\boldsymbol{\xi}} })} $

(71)

其中，系数 ${{\boldsymbol{b}}}_{}^{(n)}$ 与 ${{\boldsymbol{a}}}_{}^{(n)}$ 类似，均为n阶张量。将 $f({{\boldsymbol{\xi}} })$ 展开为Hermite多项式级数后， $f({{\boldsymbol{\xi}} })$ 关于速度 ${{\boldsymbol{\xi}} }$ 的N阶矩可表为：

$ \begin{split}& \int {f({{\boldsymbol{\xi}} }){\xi _{{i_1}}}{\xi _{{i_2}}}{\xi _{{i_3}}} \cdots {\xi _{{i_N}}}} {\mkern 1mu} {\text{d}}{{\boldsymbol{\xi}} } \\& \qquad = \int {\omega ({{\boldsymbol{\xi}} })\sum\limits_{n = 0}^N {\frac{1}{{n!}}{{\boldsymbol{a}}}_{}^{(n)} \cdot {{\boldsymbol{H}}}_{}^{(n)}({{\boldsymbol{\xi}} })} } \sum\limits_{m = 0}^N {{{\boldsymbol{b}}}_{}^{(m)} \cdot {{\boldsymbol{H}}}_{}^{(m)}({{\boldsymbol{\xi}} })} {\text{d}}{{\boldsymbol{\xi}} } { + } \\& \qquad \sum\limits_{n = N + 1}^\infty {\sum\limits_{m = 0}^N {{{\boldsymbol{a}}}_{}^{(n)}{{\boldsymbol{b}}}_{}^{(m)} \cdot \frac{1}{{n!}}} } \int {\omega ({{\boldsymbol{\xi}} }){{\boldsymbol{H}}}_{}^{(n)}({{\boldsymbol{\xi}} }){{\boldsymbol{H}}}_{}^{(m)}({{\boldsymbol{\xi}} })} {\text{d}}{{\boldsymbol{\xi}} } \end{split} $

(72)

由于Hermite多项式的正交性，式（72）右边第二项级数为0，因此有：

$ \begin{split}& \int {f({{\boldsymbol{\xi}} }){\xi _{{i_1}}}{\xi _{{i_2}}}{\xi _{{i_3}}} \cdots {\xi _{{i_N}}}} {\mkern 1mu} {\text{d}}{{\boldsymbol{\xi}} } \\& = \int {\omega ({{\boldsymbol{\xi}} })\sum\limits_{n = 0}^N {\frac{1}{{n!}}{{\boldsymbol{a}}}_{}^{(n)} \cdot {{\boldsymbol{H}}}_{}^{(n)}({{\boldsymbol{\xi}} }){p_{_N}}({{\boldsymbol{\xi}} })} } {\text{d}}{{\boldsymbol{\xi}} } = \int {{f_N}({{\boldsymbol{\xi}} }){p_{_N}}({{\boldsymbol{\xi}} }){\text{d}}{{\boldsymbol{\xi}} }} \end{split} $

(73)

其中，定义函数 $f({{\boldsymbol{\xi}} })$ 的Hermite多项式级数展开的N阶截断 ${f_N}({{\boldsymbol{\xi}} })$ （简称N阶截断式）：

$ {f_N}({{\boldsymbol{\xi}} }) \equiv \omega ({{\boldsymbol{\xi}} })\sum\limits_{n = 0}^N {\frac{1}{{n!}}{{\boldsymbol{a}}}_{}^{(n)} \cdot {{\boldsymbol{H}}}_{}^{(n)}({{\boldsymbol{\xi}} })} $

(74)

从式（73）可以看出，函数 $f({{\boldsymbol{\xi}} })$ 与其N阶截断式 ${f_N}({{\boldsymbol{\xi}} })$ 关于 ${{\boldsymbol{\xi}} }$ 的N阶矩是相等的。显然，对于 ${{\boldsymbol{\xi}} }$ 的阶数小于N的各阶矩，用 ${f_N}({{\boldsymbol{\xi}} })$ 代替 $f({{\boldsymbol{\xi}} })$ ，结果也完全一致。因此在计算 ${{\boldsymbol{\xi}} }$ 的 $n \leqslant N$ 阶矩时， ${f_N}({{\boldsymbol{\xi}} })$ 可以替代 $f({{\boldsymbol{\xi}} })$ 。

对于Boltzmann-BGK方程（式（10）），其理论上的精确解 $f(\boldsymbol{x},{{\boldsymbol{\xi}} },\;t)$ 按式（69）用Hermite多项式级数展开后，按Chapman-Enskog展开的思路， $f(\boldsymbol{x},{{\boldsymbol{\xi}} ,}\;t)$ 可以写为平衡态分布函数和一阶小量展开的和：

$ f = {f^{(0)}} + \varepsilon {f^{(1)}} + O({\varepsilon ^2}) $

(75)

$ {f^{(0)}} = \omega ({{\boldsymbol{\xi}} })\sum\limits_{n = 0}^\infty {\frac{1}{{n!}}{{\boldsymbol{a}}}_0^{(n)} \cdot {{\boldsymbol{H}}}_{}^{(n)}({{\boldsymbol{\xi}} })} $

(76)

$ {f^{(1)}} = \omega ({{\boldsymbol{\xi}} })\sum\limits_{n = 0}^\infty {\frac{1}{{n!}}{{\boldsymbol{a}}}_1^{(n)} \cdot {{\boldsymbol{H}}}_{}^{(n)}({{\boldsymbol{\xi}} })} $

(77)

注意式（76）～式（77）中的 ${{\boldsymbol{a}}}_0^{(n)}$ 与 ${{\boldsymbol{a}}}_1^{(n)}$ 为坐标排列集合的张量，下标代表Chapman-Enskog展开的 $ \varepsilon $ 阶数，而不是坐标排列。显然， $f({{\boldsymbol{x}}},{{\boldsymbol{\xi}} },\;t)$ 的各阶多项式展开系数满足：

$ {{{\boldsymbol{a}}}^{(n)}} = {{\boldsymbol{a}}}_0^{(n)} + \varepsilon {{\boldsymbol{a}}}_1^{(n)} + O({\varepsilon ^2}),n = 0,1,2, \cdots $

(78)

在此基础上，我们考察Boltzmann-BGK方程还原到各类方程需要的Hermite展开级数。

Euler方程：取 $f({{\boldsymbol{x}}},{{\boldsymbol{\xi}} },\;t) \approx {f^{(0)}}({{\boldsymbol{x}}},{{\boldsymbol{\xi}} },\;t)$ （零阶小量近似），代入式（17）～式（19），并根据式（6）的定义和 ${f^{(0)}}$ 的Maxwell-Boltzmann分布（式（11）），有：

$ {P_{ij}} = - \int {{f^{(0)}}{v_i}{v_j}{\text{d}}{{\boldsymbol{\xi}} }} = - p{\delta _{ij}} $

(79)

$ {q_i} = \frac{1}{2}\int {f{v_i}{v_j}{v_j}{\text{d}}{{\boldsymbol{\xi}} }} = \frac{1}{2}\int {{f^{(0)}}{v_i}{v_j}{v_j}{\text{d}}{{\boldsymbol{\xi}} }} = 0 $

(80)

于是式（17）～式（19）还原到Euler方程。在此过程中最高计算到 $f({{\boldsymbol{x}}},{{\boldsymbol{\xi }}},\;t)$ 的三阶矩 $ \int {f{\xi _i}{\xi _j}{\xi _j}{\text{d}}{{\boldsymbol{\xi}} }} $ （能量输运项）。根据式（73）的结论，使用三阶截断式 ${f_{N = 3}}$ 替代 $f$ ，其三阶及以下各阶矩完全相同。因此，还原到Euler方程只需要 ${f^{(0)}}$ 保留到三阶截断式。

Navier-Stokes方程：从式（79）～式（80）看到， $f$ 的零阶近似无法准确表征流体黏性应力和热传导项。因此，将 $f$ 按式（75）展开到一阶近似，同时将时间和空间导数也一并展开：

$ \frac{\partial }{{\partial t}} = \varepsilon \frac{\partial }{{\partial {t_1}}} + {\varepsilon ^2}\frac{\partial }{{\partial {t_2}}} + {O}({\varepsilon ^3}) $

(81)

$ \frac{\partial }{{\partial {x_i}}} = \varepsilon \frac{\partial }{{\partial {x_i}}} + {O}({\varepsilon ^2}) $

(82)

将分布函数 $f$ 时间和空间导数的展开式代入Boltzmann-BGK方程，有：

$ \begin{split}& \varepsilon \left[ {\frac{{\partial {f^{(0)}}}}{{\partial {t_1}}} + {{\boldsymbol{\xi}} } \cdot \nabla {f^{(0)}} + \frac{{{f^{(1)}}}}{\tau }} \right] + {\varepsilon ^2}\left[ {\frac{{\partial {f^{(0)}}}}{{\partial {t_2}}}} \right. + \\& \qquad \frac{{\partial {f^{(1)}}}}{{\partial {t_1}}} + {{\boldsymbol{\xi}} } \cdot \nabla {f^{(1)}} + \left. {\frac{{{f^{(2)}}}}{\tau }} \right] + {O}({\varepsilon ^3}) = 0 \end{split} $

(83)

忽略 ${\varepsilon ^2}$ 及更高阶量，由 ${O}(\varepsilon )$ 阶展开部分可得：

$ \frac{{{f^{(1)}}}}{\tau } + \frac{{\partial {f^{(0)}}}}{{\partial {t_1}}} + {\xi _i}\frac{{\partial {f^{(0)}}}}{{\partial {x_i}}} = 0 $

(84)

进一步分析可指出^{[37, 44]}，如果 ${f^{(1)}}$ 满足式（84）且 $ f = {f^{(0)}} + \varepsilon {f^{(1)}} $ ，则f关于 ${{\boldsymbol{\xi }}}$ 的二阶矩和三阶矩（式（4）～式（5））将准确还原Navier-Stokes方程的切应力项和热流项，进而分别还原动量和能量方程。

将 ${f^{(0)}}$ 和 ${f^{(1)}}$ 的Hermite展开式代入式（84），有：

$ \begin{split}& \omega ({{\boldsymbol{\xi}} })\sum\limits_{n = 0}^\infty {\frac{1}{{n!}}} \left[ {\frac{{{{\boldsymbol{a}}}_1^{(n)}}}{\tau } \cdot {{{\boldsymbol{H}}}^{(n)}}({{\boldsymbol{\xi}} }) + \frac{{\partial {{\boldsymbol{a}}}_0^{(n)}}}{{\partial {t_1}}} \cdot {{{\boldsymbol{H}}}^{(n)}}({{\boldsymbol{\xi}} })} \right. + \\& \qquad \left. {\frac{{\partial {{\boldsymbol{a}}}_0^{(n)}}}{{\partial {x_i}}} \cdot {\xi _i}{{{\boldsymbol{H}}}^{(n)}}({{\boldsymbol{\xi}} })} \right] = 0 \end{split} $

(85)

将式（85）的第三项所含的 ${{\xi} _i}{\boldsymbol{H}^{(n)}}({{\boldsymbol{\xi}} })$ 用递推式（62）替换，有：

$ \begin{split}& \omega ({{\boldsymbol{\xi}} })\sum\limits_{n = 0}^\infty {\frac{1}{{n!}}} \left[ {\frac{{{{\boldsymbol{a}}}_1^{(n)}}}{\tau } \cdot {{{\boldsymbol{H}}}^{(n)}}({{\boldsymbol{\xi}} }) + \frac{{\partial {{\boldsymbol{a}}}_0^{(n)}}}{{\partial {t_1}}} \cdot {{{\boldsymbol{H}}}^{(n)}}({{\boldsymbol{\xi}} })} \right. + \\& \qquad \left. {\frac{{\partial {{\boldsymbol{a}}}_0^{(n)}}}{{\partial {x_i}}} \cdot ({{\boldsymbol{H}}}_i^{(n + 1)}({{\boldsymbol{\xi}} }) + n{\delta }{{{\boldsymbol{H}}}^{(n - 1)}})} \right] = 0 \end{split} $

(86)

其中：

$ n{\delta }{{{\boldsymbol{H}}}^{(n - 1)}} = \sum\limits_{k = 1}^n {{\delta _{{i_k}i}}H_{{i_1}{i_2} \cdots {i_{k - 1}}{i_{k + 1}} \cdots {i_n}}^{(n - 1)}} $

(87)

将式（86）整理后，合并同阶Hermite多项式系数，有：

$ \begin{split}& \omega ({{\boldsymbol{\xi}} })\left[ {\sum\limits_{n = 0}^\infty {\frac{1}{{n!}}} \left( {\frac{{{{\boldsymbol{a}}}_1^{(n)}}}{\tau } + \frac{{\partial {{\boldsymbol{a}}}_0^{(n)}}}{{\partial {t_1}}}} \right)} \right. + \sum\limits_{n = 1}^\infty {\frac{1}{{n!}}} n\nabla {{\boldsymbol{a}}}_0^{(n - 1)} + \\&\qquad \left. {\quad \sum\limits_{n = - 1}^\infty {\frac{1}{{n!}}} \nabla \cdot {{\boldsymbol{a}}}_0^{(n + 1)}} \right] \cdot {{{\boldsymbol{H}}}^{(n)}}({{\boldsymbol{\xi}} }) = 0 \end{split} $

(88)

由于Hermite多项式的正交性，式（88）成立的条件为各阶 ${{{\boldsymbol{H}}}^{(n)}}({{\boldsymbol{\xi}} })$ 的系数均为0，因此有：

$ {{\boldsymbol{a}}}_1^{(n)} = - \tau \left[ {\frac{{\partial {{\boldsymbol{a}}}_0^{(n)}}}{{\partial {t_1}}} + n\nabla {{\boldsymbol{a}}}_0^{(n - 1)} + \nabla \cdot {{\boldsymbol{a}}}_0^{(n + 1)}} \right] $

(89)

式（89）指出：对于Boltzmann-BGK方程精确解的零阶和一阶近似 ${f^{(0)}}$ 和 ${f^{(1)}}$ ，其Hermite展开系数除了同阶之间满足式（78）， ${f^{(1)}}$ 的n阶展开系数 ${{\boldsymbol{a}}}_0^{(n)}$ 与 ${f^{(0)}}$ 的n+1阶展开系数 $ {{\boldsymbol{a}}}_1^{(n)} $ 之间还存在依赖关系。我们在 ${f^{(0)}}$ 的 ${N_0}$ 阶Hermite截断式 $ {f}_{{N}_{0}}^{(0)} $ 的基础上构造Boltzmann-BGK方程的近似解 $\tilde f$ ：

$ \tilde{f} = {f}_{{N}_{0}}^{(0)}+\epsilon {\tilde{f}}^{(1)} = {\tilde{f}}^{(0)}+\epsilon {\tilde{f}}^{(1)} $

(90)

其中， ${\tilde f^{(0)}} = f_{{N_0}}^{(0)}，$ ${N_0}$ 将在以下的分析中确定。将 ${\tilde f^{(0)}}$ 和 ${\tilde f^{(1)}}$ 用Hermite多项式级数展开：

$ {\tilde f^{(0)}} = \omega \sum\limits_{n = 0}^\infty {\frac{1}{{n!}}{\tilde {\boldsymbol{a}}}_0^{(n)} \cdot {{{\boldsymbol{H}}}^{(n)}}({{\boldsymbol{\xi}} })} $

(91)

$ {\tilde f^{(1)}} = \omega \sum\limits_{n = 0}^\infty {\frac{1}{{n!}}{\tilde {\boldsymbol{a}}}_1^{(n)} \cdot {{{\boldsymbol{H}}}^{(n)}}({{\boldsymbol{\xi}} })} $

(92)

显然 $n \leqslant {N}_{0}$ 时 ${\tilde {\boldsymbol{a}}}_0^{(n)} = {{\boldsymbol{a}}}_0^{(n)}$ ； $n > {N_0}$ 时 ${\tilde {\boldsymbol{a}}}_0^{(n)} = 0$ 。为了方便后续描述，这里依然保留 ${\tilde f^{(0)}}$ 和 ${\tilde {\boldsymbol{a}}}_0^{(n)}$ 标记。可以证明， ${\tilde {\boldsymbol{a}}}_1^{(n)}$ 和 ${\tilde {\boldsymbol{a}}}_0^{(n)}$ 的关系类似于 ${{\boldsymbol{a}}}_1^{(n)}$ 和 ${{\boldsymbol{a}}}_0^{(n)}$ （附录A中式（131）和式（132））。如需 ${\tilde f^{(1)}}$ 对 ${{\boldsymbol{\xi}} }$ 的 $N$ 阶矩与精确解 ${f^{(1)}}$ 对 ${{\boldsymbol{\xi}} }$ 的 $N$ 阶矩相同，则需要 $\tilde f_N^{(1)} = f_N^{(1)},$ 即 ${\tilde f^{(1)}}$ 与 ${f^{(1)}}$ 的 $N$ 阶截断式相等（但这不意味着 ${\tilde f^{(1)}} = f_N^{(1)}$ ），因此 $\tilde{\boldsymbol{a}}_{1}^{(n)} = {{{\boldsymbol{a}}}}_{1}^{(n)}, 0 \leqslant n\leqslant N $ 。此时根据式（132）约束， ${\tilde {\boldsymbol{a}}}_1^{(N)}$ 由 ${{\boldsymbol{a}}}_0^{(N - 1)}$ ～ ${{\boldsymbol{a}}}_0^{(N + 1)}$ 确定，因而要求

$\tilde{\boldsymbol{a}}_{0}^{(n)} = {{{\boldsymbol{a}}}}_{0}^{(n)},\;\;\;\;0\leqslant n \leqslant N+1 $

(93)

即 $ {\tilde{f}}^{(0)} = {f}_{N+1}^{(0)} $ （或 $ {N_0} = N + 1 $ ）， ${\tilde f^{(0)}}$ 为平衡态分布函数 ${f^{(0)}}$ 的 $N + 1$ 阶截断式。

根据上述讨论，在Chapman-Enskog展开一阶小量近似下，如需 $\tilde f$ 对 ${\boldsymbol\xi }$ 的 $N$ 阶矩与f对 ${{\boldsymbol{\xi}} }$ 的 $N$ 阶矩相等，要求 ${\tilde f^{(0)}}$ 、 ${\tilde f^{(1)}}$ 对 ${{\boldsymbol{\xi}} }$ 的 $N$ 阶矩分别与精确解 ${f^{(0)}}$ 、 ${f^{(1)}}$ 对 ${{\boldsymbol{\xi}} }$ 的 $N$ 阶矩相同；对于 ${\tilde f^{(0)}}$ 而言，要求 ${\tilde f^{(0)}} = f_{N + 1}^{(0)}$ ，进而 $ {\tilde{\boldsymbol{ a}}}_0^{(N + 1)} = {{\boldsymbol{a}}}_0^{(N + 1)} \ne 0 $ ，而式（133）则将 $ {\tilde{\boldsymbol{ a}}}_0^{(N + 1)} $ 与 $ {\tilde {\boldsymbol{a}}}_1^{(N + 2)} $ 联系起来，这意味着对于 ${\tilde f^{(1)}}$ ，至少有：

$ {\tilde {\boldsymbol{a}}}_1^{(N + 2)} \ne 0 $

(94)

这也就是说， ${\tilde f^{(1)}}$ 的多项式展开级数至少为N+2阶，尽管 $ {\tilde {\boldsymbol{a}}}_1^{(N + 2)} \ne {{\boldsymbol{a}}}_1^{(N + 2)} $ ，但 $ {\tilde {\boldsymbol{a}}}_1^{(N + 2)} $ 的非零性使得 ${\tilde f^{(1)}}$ 至少为 $N + 2$ 阶Hermite多项式级数（但 ${\tilde f^{(1)}}$ 不是精确解 ${f^{(1)}}$ 的 $N + 2$ 阶截断式 $f_{N + 2}^{(1)}$ ），因而近似解 $\tilde f$ 也至少为 $N + 2$ 阶Hermite多项式级数。

在以上讨论的基础上，我们可以回到具体需还原的方程类型中考察N的值。对于等温弱可压Navier-Stokes方程，只需要准确还原质量和动量方程，故需 $\tilde f$ 对 ${{\boldsymbol{\xi}} }$ 的二阶矩与f对 ${{\boldsymbol{\xi}} }$ 的二阶矩相同，要求 $ {\tilde{f}}_{2}^{(1)} = {f}_{2}^{(1)} $ ，于是N = 2，根据式（93）有 ${\tilde f^{(0)}} = f_3^{(0)}$ ；式（94）显示， ${\tilde f^{(1)}}$ 实际为四阶Hermite多项式级数，因此 $\tilde f$ 也至少为四阶Hermite多项式级数；于是计算 $\tilde f$ 对 ${{\boldsymbol{\xi}} }$ 的二阶矩时，积分内核 $ \tilde f{\xi _i}{\xi _j}/\omega $ 至少为六阶多项式。类似的，对于可压缩流动的Navier-Stokes方程，需要能量方程也能准确还原，因此需要 $\tilde f$ 对 ${{\boldsymbol{\xi}} }$ 的三阶矩与f对 ${{\boldsymbol{\xi}} }$ 的三阶矩相同，故N = 3， ${\tilde f^{(0)}} = f_4^{(0)}$ ；根据式（94）， ${\tilde {\boldsymbol{a}}}_1^{(5)} \ne 0$ ，因此 ${\tilde f^{(1)}}$ 和 $\tilde f$ 也至少为五阶Hermite多项式级数， $ \widetilde f{\xi _i}{\xi _j}{\xi _j}/\omega $ 至少为八阶多项式，这与Shan等^[30]的结论一致。

关于积分内核的多项式阶数，也可以从另一个方面去解释^[45]：对于等温弱可压Navier-Stokes方程，需 $\tilde f$ 对 ${{\boldsymbol{\xi}} }$ 的二阶矩与f对 ${{\boldsymbol{\xi}} }$ 的二阶矩相同，因此要求 $ {\tilde{f}}_{2}^{(1)} = {f}_{2}^{(1)} $ ，于是N = 2，根据式（93）有 ${\tilde f^{(0)}} = f_3^{(0)}$ 。要想在数值计算过程中准确表达 $f_3^{(0)}$ ，等价于需准确表达 $f_3^{(0)}$ 的三阶Hermite展开式系数 $ {{\boldsymbol{a}}}_0^{(3)} $ ，这需要计算 $f_3^{(0)}$ 对 ${{\boldsymbol{\xi}} }$ 的三阶矩，积分内核 $ f_3^{(0)}{\xi _i}{\xi _j}{\xi _k}/\omega $ 为六阶多项式；类似的，对于可压缩流动的Navier-Stokes方程，有 ${\tilde f^{(0)}} = f_4^{(0)}$ ，要想在数值计算过程中准确表达 $f_4^{(0)}$ ，等价于需准确表达 $f_4^{(0)}$ 的四阶Hermite展开式系数 ${{\boldsymbol{a}}}_0^{(4)}$ ，此时需计算 $f_4^{(0)}$ 对 ${{\boldsymbol{\xi}} }$ 的四阶矩，因此积分内核 $ f_4^{(0)}{\xi _i}{\xi _j}{\xi _k}{\xi _l}/\omega $ 为八阶多项式。

综上，我们将具体还原的方程类型和其多项式阶数列在表1中。

${f^{(0)}}$ 的零到四阶Hermite展开式系数为：

$ {{\boldsymbol{a}}}_0^{(0)} = \rho ,\quad {{\boldsymbol{a}}}_0^{(1)} = \rho {{\boldsymbol{u}}} $

(95)

$ {{\boldsymbol{a}}}_0^{(2)} = \rho\big [{{{\boldsymbol{u}}}^2} + (\theta - 1)\boldsymbol{\delta} \big] $

(96)

$ {{\boldsymbol{a}}}_0^{(3)} = \rho \big[{{{\boldsymbol{u}}}^3} + (\theta - 1){\boldsymbol{\delta} {\boldsymbol{u}}}\big] $

(97)

$ {{\boldsymbol{a}}}_0^{(4)} = \rho \big[{{{\boldsymbol{u}}}^4} + (\theta - 1){{\boldsymbol{\delta}} }{{{\boldsymbol{u}}}^2} + {(\theta - 1)^2}{{{\boldsymbol{\delta}} }^2}\big] $

(98)

其中，矢量或张量并列（ $ {{{\boldsymbol{u}}}^2} $ 、 $ {{{\boldsymbol{u}}}^3} $ 、 $ {{{\boldsymbol{\delta}} }^2} $ 等）或者矢量与张量并列（ $ {{\boldsymbol{\delta}} {\boldsymbol{u}}} $ 、 $ {{\boldsymbol{\delta}} }{{{\boldsymbol{u}}}^2} $ 等），为张量对称量表示法^[43]。 ${f^{(0)}}$ 的二到四阶Hermite截断式可表达为：

$\qquad\qquad \begin{split} f_2^{(0)}({{\boldsymbol{\xi}} }) =& \rho \omega ({{\boldsymbol{\xi}} })\Big\{ 1 + {{\boldsymbol{\xi}} } \cdot {{\boldsymbol{u}}}+ \\& \frac{1}{2}\big[{({{\boldsymbol{\xi}} } \cdot {{\boldsymbol{u}}})^2} - {u^2} + (\theta - 1)({\xi ^2} - D)\big]\Big\} \end{split} $

(99)

$\qquad\qquad \begin{split} f_3^{(0)}({{\boldsymbol{\xi}} }) =& f_{2}^{(0)}({{\boldsymbol{\xi}} }) + \rho \omega ({{\boldsymbol{\xi}} })\bigg\{ \frac{{{{\boldsymbol{\xi}} } \cdot {{\boldsymbol{u}}}}}{6}\Big[{({{\boldsymbol{\xi}} } \cdot {{\boldsymbol{u}}})^2} - \\&3{u^2} + 3(\theta - 1)({\xi ^2} - D - 2)\Big]\bigg\} \end{split} $

(100)

$ \begin{split}& f_4^{(0)}({{\boldsymbol{\xi}} }) = f_{3}^{(0)}({{\boldsymbol{\xi}} }) + \rho \omega (\boldsymbol{\xi })\cdot \\& \Bigg\{\frac{1}{{24}}\Big[{(\boldsymbol{\xi } \cdot {{\boldsymbol{u}}})^4} - 6{(\boldsymbol{\xi } \cdot {{\boldsymbol{u}}})^2}{u^2} + 3{u^4}\Big] + \\& \frac{{\theta - 1}}{4}\Big\{ ({\xi ^2} - D - 2)\Big[{(\boldsymbol{\xi } \cdot {{\boldsymbol{u}}})^2} - {u^2}\Big] - 2{(\boldsymbol{\xi } \cdot {{\boldsymbol{u}}})^2}\Big\} +\\& \frac{{{{(\theta - 1)}^2}}}{8}\Big[{\xi ^4} - 2(D + 2){\xi ^2} + D(D + 2)\Big]\;\Bigg\} \end{split} $

(101)

在表1的指导下，可根据实际计算的流动类型选用对应阶数的 ${f^{(0)}}$ 截断式。下一步还需要将 $\tilde f$ 和 ${\tilde f^{(0)}}$ 的截断式以及它们对 ${{\boldsymbol{\xi}} }$ 的各阶矩离散化，即使用合理的数值方法精确计算积分。为了不失一般性，对于具有K阶精度的数值积分公式可写为：

$ \int {\omega ({{\boldsymbol{\xi}} })g({{\boldsymbol{\xi}} })} {\mkern 1mu} {\text{d}}{{\boldsymbol{\xi}} = }\sum\limits_{a = 1}^d {{w_a}g(} {{{\boldsymbol{\xi}} }_a}) $

(102)

这里K阶精度的意义是对于 $ g({{\boldsymbol{\xi}} }) $ 为小于或等于K阶多项式，式（102）都能精确成立。该式中的积分点 $ {{{\boldsymbol{\xi}} }_a} $ 与权重 $ {w_a} $ 如何选取将在稍后讨论。

定义分布函数的离散形式：

$ {f_a} = \frac{{{w_a}\tilde f({{{\boldsymbol{\xi}} }_{{\boldsymbol{a}}}})}}{{\omega \left( {{{{\boldsymbol{\xi}} }_{{\boldsymbol{a}}}}} \right)}},\;f_a^{(0)} = \frac{{{w_a}f_{N + 1}^{(0)}({{{\boldsymbol{\xi}} }_{{\boldsymbol{a}}}})}}{{\omega \left( {{{{\boldsymbol{\xi}} }_{{\boldsymbol{a}}}}} \right)}} $

(103)

在Boltzmann-BGK方程（10）中用近似解 $ \tilde f $ 、 $ f_{N + 1}^{(0)} $ 分别替换f与 $ f_{}^{(0)} $ ，代入 ${{\boldsymbol{\xi}} = }{{{\boldsymbol{\xi }}}_a}$ 并在方程两边同时乘以权重 $ {w_a} $ 后，LBM方程可简洁表示为：

$ \frac{{\partial {f_a}}}{{\partial t}} + {{{\boldsymbol{\xi}} }_{a}} \cdot \nabla {f_a} = - \frac{{{f_a} - f_a^{(0)}}}{\tau } $

(104)

附录B中我们将证明，只要 $ f_{}^{(0)} $ 的截断式阶数满足表1，且 $ K\geqslant M+N $ ，则式（104）可还原到对应的Euler方程、等温弱可压Navier-Stokes方程或可压缩热流Navier-Stokes方程。

在离散化 $ f_{N + 1}^{(0)} $ 时，根据之前的论述，满足表1 的 $ f_{N + 1}^{(0)} $ 将使得近似解 $ \tilde f $ 对 ${{\boldsymbol{\xi}} }$ 的各阶矩与精确解的计算结果一致。根据式（1）～式（2）和式（4），以及积分公式（102），密度、速度和温度可用离散分布函数表示为：

$ \begin{split} \rho =& \int {f\;{\text{d}}{{\boldsymbol{\xi}} }} = \int {\tilde f\;{\text{d}}{{\boldsymbol{\xi}} }} = \int {\omega (\xi )\frac{{\tilde f\;}}{{\omega (\xi )}}{\text{d}}{{\boldsymbol{\xi}} }} \\ = &\sum\limits_{a = 1}^d {\frac{{{w_a}\tilde f({\xi _a})\;}}{{\omega ({\xi _a})}}} = \sum\limits_{a = 1}^d {{f_a}} \end{split} $

(105)

$ \rho u = \int {\tilde f{{\boldsymbol{\xi}} }\;{\text{d}}{{\boldsymbol{\xi}} }} = \sum\limits_{a = 1}^d {\frac{{{w_a}\tilde f({\xi _a}){{{\boldsymbol{\xi}} }_a}\;}}{{\omega ({\xi _a})}}} = \sum\limits_{a = 1}^d {{f_a}{{{\boldsymbol{\xi}} }_a}} $

(106)

$ \begin{split} \rho {u^2} + pD = &\rho ({u^2} + \theta D) = \int {\tilde f{\xi _i}{\xi _i}{\text{d}}{{\boldsymbol{\xi}} }} \\=& \sum\limits_{a = 1}^d {\frac{{{w_a}\tilde f({\xi _a})\xi _a^2\;}}{{\omega ({\xi _a})}}} = \sum\limits_{a = 1}^d {{f_a}\xi _a^2} \end{split} $

(107)

在等温弱可压流动中， $\theta = 1$ ；根据式（21）和式（61）关于 $ {w_a} $ 和 $ \omega (\boldsymbol\xi ) $ 的定义，有 $ {w_a} = \omega ({\boldsymbol{\xi }_a}) $ ，于是 ${f_a} = \tilde f({\xi _a})，$ 这使得式（105）～式（106）与式（40）～式（41）一致。同时，在式（103）中将 $ {\boldsymbol{\xi }_a} $ 与权重 $ {w_a} $ 代入 $ f_{N + 1}^{(0)} $ 后，与等温弱可压流动类似，可将 $ f_{N + 1}^{(0)} $ 中的 $ {{\boldsymbol{\xi}}}_{a} $ 、 ${{\boldsymbol{u}}}$ 、 $\theta $ 分别利用关系式 ${{{\boldsymbol{c}}}_a} = {c_s}{\boldsymbol{\xi }_a}$ 、 ${{\boldsymbol{u}}}' = {c_s}{{\boldsymbol{u}}}$ 、 $\theta ' = c_s^2\theta $ 替换为归一化的 ${{{\boldsymbol{c}}}_a} $ 、 ${{\boldsymbol{u}}}' $ 、 $\theta '$ ，以方便计算程序的编制；其中，归一化常数 ${c_s} = 1/c$ 将在离散速度模型的求解中确定。

4.2 时间和空间离散

多层速度模型的离散Boltzmann方程（式（104））与等温弱可压流动相同，因此可以沿用2.2节类似的思路，在特征线 ${{{\boldsymbol{\xi}} }_a}$ 上分别积分（式（28）），并使用二阶时间格式（式（31））对时间积分离散化。此时与等温流动的不同之处在于，流体黏性 $\mu $ 成为了与温度变化相关的量，例如气体的粘性随温度的变化服从Sutherland公式：

$ \frac{\mu }{{{\mu _0}}} = {\left( {T/{T_0}} \right)^{3/2}}\frac{{1 + S/{T_0}}}{{T/{T_0} + S/{T_0}}} $

(108)

其中，S为随气体类型变化的常数，对于空气S = 110.4 K； $ {\mu _0} $ 为气体在温度 $ {T_0} $ 时的动力学黏性。此时松弛时间 $ \bar \tau $ 不再是定值，而是与压强（温度）温度相关的变量：

$ \bar \tau - \frac{{\Delta t}}{2} = \frac{\mu }{p} $

(109)

对于给定雷诺数Re和马赫数Ma的可压缩热流问题（以及指定了定义雷诺数、马赫数的尺度L和速度U），如果在尺度 $L$ 上分布有 ${N_L}$ 个均匀LB网格，则式（12）中的参考量 ${L_0}$ 为 ${L_0} = L/{N_L}$ ，对应的无量纲参考松弛时间为：

$ \bar \tau _0^{} = \frac{{{N_L}{{Ma}}}}{{{c_s}{{Re}} }} + \frac{1}{2} $

(110)

事实上式（110）也适用于等温弱可压流动。对于流场中任意一点的分布函数 ${f_a}({{\boldsymbol{x}}},t)$ ，先根据式（105）～式（107）计算出 $p$ 和 $\rho $ ，该点的温度为：

$ \frac{T}{{{T_0}}} = \frac{{p/\rho }}{{p_0^{}/\rho _0^{}}} = p/\rho $

(111)

代入式（108）中求得 $ \mu /{\mu _0} $ ，再将 $ \mu /{\mu _0} $ 、 $p/{p_0} = p$ 代入式（109），得到：

$ \bar \tau = \frac{\mu }{{{\mu _0}p}}\left( {\bar \tau _0^{} - \frac{1}{2}} \right) + \frac{1}{2} $

(112)

即为适应于温度变化的松弛时间。

4.3 离散速度模型

在之前的讨论中，有一个重要的问题没有展开——如何获得具有K阶精度的数值积分公式（式（102））。我们指出，对于任意的m阶多项式 $g(\xi ) = {p_m}(\xi )$ ，数值积分公式（102）成立的充要条件是：

$ \sum\limits_{a = 1}^d {{w_a}H_{j}^{(i)}(} {\xi _a}) = {\delta _{0i}},\quad i = 0,1,2, \cdots ,m $

(113)

其证明见附录C，求解该式即可获得 $ {\boldsymbol{\xi }_a} $ 和 ${w_a}$ 。特别说明：在多维情况下，下标排列 ${{\boldsymbol{j}}} = ({j_1},{j_2},{j_3}, \cdots ,{j_i})$ 是长度为i的任意排列，此时i阶多项式 $ H_{j}^{(i)} $ 有多个，因此对于给定的i，式（113）是多个方程。例如二维情况下，i = 2时，式（113）实际为：

$ \begin{array}{l}{\displaystyle \sum _{a = 1}^{d}{w}_{a}{H}_{11}^{(2)}(}{\xi }_{a}) = 0,\;\;\;\;{\displaystyle \sum _{a = 1}^{d}{w}_{a}{H}_{12}^{(2)}(}{\xi }_{a}) = 0,\\ {\displaystyle \sum _{a = 1}^{d}{w}_{a}{H}_{22}^{(2)}(}{\xi }_{a}) = 0\end{array} $

为了求解式（113），我们预设解的四个性质：

（1）均匀性。一维情况下， ${\xi}_a$ 序列总是某一特征常数c （ ${c_s} = 1/c$ ）的整数倍，即 $\{ {\xi _a}\} $ = $\{ 0, \pm c,$ $ \pm 2c, \cdots \pm nc\} ,$ n为待定常数；多维情况下 $ {{\boldsymbol\xi }_a} $ 的每个方向分量都是c的整数倍。由于LBM一般采用均匀网格，这样可以使得沿 $ {\boldsymbol{\xi }_a} $ 方向特征线积分 $\Delta t$ 后， $\boldsymbol{x}+ \boldsymbol{\xi }_a\Delta t$ （即归一化的 $\boldsymbol{x}+\boldsymbol{c}_a\Delta t$ ）依然落在网格点上（on-lattice），避免插值等操作带来额外误差。

（2）对称性。关于坐标轴对称和原点对称的 $ {\boldsymbol{\xi }_a} $ 对应的权重 ${w_a}$ 相同，多维情况下 $ {\boldsymbol{\xi }_a} $ 各方向分量重新排列对应的权重 ${w_a}$ 也相同。具体而言，对于整数i、j和k，一维情况下有 ${w_a}( \pm ic) = {w_i}$ ；二维情况下有 ${w_a}( \pm ic, \pm jc) = {w_a}( \pm jc, \pm ic) = {w_{ij}}$ ；三维情况下有 $ {w}_{a}(\pm ic,\pm jc\pm kc) = {w}_{a}(\pm jc,\pm ic,\pm kc) $ $ = \cdots = $ ${w}_{a}(\pm kc，\pm jc,\pm ic) = {w}_{ijk}，$ 共6种 $i、j、k$ 排列顺序。

（3）正定性。权重 ${w_a}$ 需为非负值，才可以保证计算的稳定性。

（4）经济性。对于同样阶数的积分精度可能存在多个解，我们倾向于选择积分点数d最少的解，这样可以最大化节约计算资源。

在这些预设条件下，求解式（113）过程分为三步：

（1）根据积分精度确定方程。奇数阶的Hermite多项式至少是 $ \boldsymbol{\xi } $ 某一个坐标分量的奇函数，由于积分点和权重的对称性，奇数阶Hermite多项式是无条件满足式（113）的，对于求解 $ {\boldsymbol{\xi }_a} $ 和 ${w_a}$ 没有贡献，故只有偶数阶Hermite多项式才能贡献有效的方程；并且，如果 $ {\boldsymbol{\xi }_a} $ 和 ${w_a}$ 使得式（113）满足2k （k>0）阶Hermite多项式，则该式对2k+1阶Hermite多项式也自动成立，即积分精度实际具有2k+1阶。因此，要求解2k+1阶精度的积分公式，先统计 $ {{H}^{(0)}}({\boldsymbol\xi }) $ ～ $ {{H}^{(2k)}}(\boldsymbol{\xi }) $ 中偶数阶、不具备重复对称性的多项式个数p，则式（113）可列出p个方程，可求解p个未知数（权重）。按经济性要求，我们一般按 ${\boldsymbol{\xi }_a}$ 的模大小 $|{{\boldsymbol\xi }_a}| = {({\boldsymbol{\xi }_a} \cdot {\boldsymbol{\xi }_a})^{1/2}}$ 对离散点进行分组，选择 ${\boldsymbol{\xi }_a}$ 的离散点数量最少、且 $|{\boldsymbol{\xi }_a}|$ 最小的前p组对应的 ${w_a}$ 作为这些方程的未知数。但这种方案可能导致方程组无解，因此需要剔除前p组的部分未知数，加入 ${\boldsymbol{\xi }_a}$ 数量次少的离散点组构成方程。此过程一般需要穷举法逐一验证。

（2）求解步骤（1）中的p个方程，得到 ${w_a}$ 与c的关系。

（3）确定特征常数c。只要所有 ${w_a}$ 非负，c可以取任意正实数，我们可以选择适当的c，使得步骤（2）解出的一个或者多个 ${w_a}$ 为零，其余的 ${w_a}$ 均为正值，进而得到离散点最少的积分点公式。

该求解过程已在Shan^[30]中初步给出了思路，并由Nie^[45]和Shan^[46-47]对求解过程进行了更详细地描述。Spiller等^[48-50]则将求解方法编制成了Python程序包。以下在一维和二维情况下以实例说明该过程。

一维情况下，我们考察一到九阶精度的积分公式。对求解 ${\xi _a}$ 和 ${w_a}$ 有贡献的Hermite多项式为 ${H}^{(0)}(\xi )$ 、 ${H}^{(2)}(\xi )$ 、··· ${H^{(8)}}(\xi )$ ，可列出5个方程。同时，选择未知数的个数不少于方程数量，例如 ${w_0}～{w_5}$ ，共6个未知权重，对应积分点 $\left|{\boldsymbol\xi }_{a}\right|\leqslant5c$ 。将有贡献的多项式、选取的积分点和权重代入式（113），并对各方程依次高斯消元（可借助Mathematica、Sympy等计算机代数系统完成），得：

$ \begin{array}{l} \begin{array}{ccc cccc} \quad \quad \quad & {w}_{0}& {w}_{1}& {w}_{2}& {w}_{3}& {w}_{4}& {w}_{5} \end{array}\\ \begin{array}{c} {H}^{(0)}:\\ {H}^{(2)}:\\ {H}^{(4)}:\\ {H}^{(6)}:\\ {H}^{(8)}: \end{array} \left[ \begin{array}{cccccc} 1& ~~2& ~~2& ~~2& ~~2& 2\\ 0& ~~1& ~~4& ~~9& ~~16& 25\\ 0& ~~0& ~~1& ~~6& ~~20& 50\\ 0& ~~0& ~0& ~~1& ~~8& 35\\ 0& ~~0& ~0& ~~0& ~~1& 10 \end{array} \right] \cdot {{\boldsymbol{w}}} = {{\boldsymbol{r}}}\\ \end{array} $

(114)

其中， ${w}= {({w_0},{w_1},{w_2}, \cdots ,{w_5})^{\text{T}}}$ 。

右端项 $ \boldsymbol{r} = $ $({r_1},{r_2},$ ${r_3},{r_4}, {r_5}{)^{\text{T}}}$ 为：

$\boldsymbol {r} = \left[ \begin{array}{l} 1 \\ 1/2{c^2} \\ - ({c^2} - 3)/24{c^2} \\ (4{c^4} - 15{c^2} + 15)/720{c^6} \\ - (12{c^6} - 49{c^4} + 70{c^2} - 35)/13440{c^8} \\ \end{array} \right] $

(115)

为了方便对照，我们在式（114）系数矩阵的左侧和顶端分别列出了各系数对应的Hermite多项式和未知数。如需要五阶积分公式，需要考察 ${H^{(0)}}(\xi )～$ ${H^{(4)}}(\xi )$ 对应的前3个方程。按照一般情况下未知数与方程数相当的原则，取未知数为 ${w_0}～{w_2}，$ 对应的积分点为 $ {\xi }_{a} = \{0，\pm c，\pm 2c\}，$ 其余的 ${w_3}～{w_5}$ 均设置为0（从这里可以看到式（114）的未知数只要不少于有效个数，其数量选择具有一定的任意性）。求解这3个方程得：

$\begin{split} & {w_0} = (7{c^2} - {c^4} - 4)/(4{c^2})\\ & {w_1} = ({c^4} - 3{c^2} + 3)/(6{c^2})\\ & {w_2} = (3 - {c^2})/{24} \end{split} $

(116)

此时，只要给定任意使 ${w_a}$ 满足正定性的 $c > 0$ ，按式（116）的 ${w_0}～{w_2}$ 即可计算五阶精度的积分。如果取 $c = \sqrt 3 $ ，则可以使得 ${w_2} = 0$ ，进而使用三个积分点 $ {\xi }_{a} = \{0，\pm \sqrt{3}\} $ 也能获得具有五阶精度（最经济）的积分公式。这与使用马赫展开法求解得到的D1Q3模型是一致的。

进一步地，我们考虑使用积分点 $ \left|{\xi }_{a}\right|\leqslant 2c $ 能否计算七阶积分。此时需要式（114）的前4个方程。令 ${w_3}～{w_5}$ 均为0，则 ${H^{(6)}}$ 对应的方程将成为 $0 = {r_4}$ ，但 ${r_4}$ 是恒为正的，因而方程无解。这说明使用积分点 $ \left|{\xi }_{a}\right|\leqslant2c $ 不可能得到有意义的七阶积分公式。增加未知数到 ${w_0}～{w_3}$ （ ${w_4} = {w_5} = 0$ ），相应的积分点为 $ {\xi }_{a} = \{0，\pm c，\pm 2c，\pm 3c\}，$ 可解得：

$ \begin{split}& {w_0} = (36{c^6} - 49{c^4} + 42{c^2} - 15)/36{c^6} \\& {w_1} = (12{c^4} - 13{c^2} + 5)/16{c^6} \\& {w_2} = ( - 3{c^4} - 10{c^2} + 5)/40{c^6} \\& {w_3} = (4{c^4} - 15{c^2} + 15)/720{c^6} \end{split} $

(117)

此时 ${w_1}$ 与 ${w_3}$ 都是恒正的，为了得到最少积分点，只能取 ${w_2} = 0$ ，得到 ${c^2} = (5 \pm \sqrt {10} )/3$ ，代入式（117）得到 ${w_0}$ 、 ${w_1}$ 和 ${w_3}$ 的具体数值，对应的积分点组为 $ {\xi }_{a} = \{0，\pm c，\pm 3c\} $ 。进一步考虑 $ \left|{\xi }_{a}\right|\leqslant 3c $ 能否计算九阶积分，需要 ${H^{(0)}}～{H^{(8)}}$ 对应的5个方程。此时由于 $ {w}_{4} = {w}_{5} = 0，$ ${H^{(8)}}$ 对应的方程为 $ 0 = {r}_{5}，$ 幸运的是该约束条件可解出 ${c^2} \approx 1.19698$ （否则需要考虑 ${w_4} \ne 0$ 的情况）。由于前4个方程和未知数与七阶积分情况没有变化，可代入式（117）得到对应的 ${w_0}～{w_3}$ 。以此类推，可得到更高阶数的最经济积分点与权重^[47]，此处不再赘述。

多维情况下，对于2n阶Hermite多项式 $ {H}_{\boldsymbol{i}}^{(2n)} $ ，由于坐标分量的对称性，对求解积分公式有贡献的多项式，只需要考察排列：

$\boldsymbol{i} = (\underset{2\alpha 个}{\underbrace{1,\cdots 1}},\underset{2\beta 个}{\underbrace{2,\cdots 2}},\underset{2\gamma 个}{\underbrace{3,\cdots 3}}) = ({1}^{2\alpha }{2}^{2\beta }{3}^{2\gamma }) $

(118)

其中， $ \alpha $ 、 $\;\beta $ 、 $\gamma $ 均为整数，满足 $2\alpha + 2\beta + 2\gamma = 2n$ ，且 $\alpha \geqslant \beta \geqslant \gamma \geqslant0，$ 最右端为本文约定的简写形式。

二维情况下，对于一到九阶精度积分公式，有效的Hermite多项式依次为 $ {H}^{(0)} $ 、 ${H}_{11}^{(2)}$ 、 $ {H}_{{1}^{4}}^{(4)} $ 、 $ {H}_{{1}^{2}{2}^{2}}^{(4)} $ 、 $ {H}_{{1}^{6}}^{(6)} $ 、 $ {H}_{{1}^{4}{2}^{2}}^{(6)} $ 、 $ {H}_{{1}^{8}}^{(8)} $ 、 $ {H}_{{1}^{6}{2}^{2}}^{(8)} $ 、 $H_{{1^4}{2^4}}^{(8)}$ 。由一维情况的结果启发，我们尝试用 $ \left|{\xi }_{ai}\right|\leqslant3c $ 的积分点组和权重来构造积分公式，如表2所示。将有效Hermite多项式和积分点、权重代入式（113），得到与式（114）类似、经过高斯消元的方程组 $A \cdot {{\boldsymbol{w}}} = {{\boldsymbol{r}}},$ 其中系数矩阵A为：

$ \begin{array}{l} \begin{array}{ccccc ccccc c} \quad\quad\quad\quad& {w}_{0}& {w}_{10}& {w}_{11}& {w}_{20}& {w}_{21}& {w}_{22}& {w}_{30}& {w}_{31}& {w}_{32}& {w}_{33} \end{array}\\ \begin{array}{c} {H}^{(0)}:\\ {H}_{11}^{(2)}:\\ {H}_{{1}^{2}{2}^{2}}^{(4)} :\\ {H}_{{1}^{4}}^{(4)}:\\ {H}_{{1}^{4}{2}^{2}}^{(6)} :\\ {H}_{{1}^{6}}^{(6)}:\\ {H}_{{1}^{4}{2}^{4}}^{(8)} :\\ {H}_{{1}^{6}{2}^{2}}^{(8)} :\\ {H}_{{1}^{8}}^{(8)}: \end{array} \left[ \begin{array}{ccccc ccccc} ~~~1& ~~~~4& ~~~~4& ~~~~4& ~~~8& ~~4& ~~~4& ~~~8& ~~8& ~4\\ ~~~0& ~~~~1& ~~~~2& ~~~~4& ~~~10& ~~8& ~~~9& ~~~20& ~~26& ~18\\ ~~~0& ~~~~0& ~~~~1& ~~~~0& ~~~8& ~~16& ~~~0& ~~~18& ~~72& ~81\\ ~~~0& ~~~~0& ~~~~0& ~~~~1& ~~~2& ~~2& ~~~6& ~~~12& ~~14& ~12\\ ~~~0& ~~~~0& ~~~~0& ~~~~0& ~~~1& ~~4& ~~~0& ~~~6& ~~33& ~54\\ ~~~0& ~~~~0& ~~~~0& ~~~~0& ~~~0& ~~0& ~~~1& ~~~2& ~~2& ~2\\ ~~~0& ~~~~0& ~~~~0& ~~~~0& ~~~0& ~~1& ~~~0& ~~~0& ~~12& ~36\\ ~~~0& ~~~~0& ~~~~0& ~~~~0& ~~~0& ~~0& ~~~0& ~~~1& ~~4& ~9\\ ~~~0& ~~~~0& ~~~~0& ~~~~0& ~~~0& ~~0& ~~~0& ~~~0& ~~0& ~0 \end{array} \right] \end{array} $

右端项r为：

$ \boldsymbol{r} = \left[ \begin{array}{l} 1 \\ 1/2{c^2} \\ 1/4{c^4} \\ (3 - {c^2})/24{c^4} \\ (3 - {c^2})/48{c^6} \\ (4{c^4} - 15{c^2} + 15)/720{c^6} \\ {({c^2} - 3)^2}/576{c^8} \\ (4{c^4} - 15{c^2} + 15)/1440{c^8} \\ 12{c^6} - 49{c^4} + 70{c^2} - 35 \end{array} \right] $

(119)

对于五阶积分公式，需满足 $ {H}^{(0)} $ 、 ${H}_{11}^{(2)}$ 、 $ {H}_{{1}^{4}}^{(4)} $ 、 $H_{{1^2}{2^2}}^{(4)}$ 对应的前4个方程。我们按数量最少原则，取前4个积分点组的未知数 $ {w}_{0} $ 、 $ {w}_{10} $ 、 $ {w}_{11} $ 、 $ {w}_{20} $ ，其余未知数均设置为0，求解得到：

$\begin{split}& {w_0} = 1 - \frac{5}{{2{c^2}}} + \frac{5}{{2{c^4}}} \\& {w_{10}} = \frac{{2{c^2} - 3}}{{3{c^4}}} \\& {w_{11}} = \frac{1}{{4{c^4}}} \\& {w_{20}} = - \frac{{{c^2} - 3}}{{24{c^4}}} \end{split} $

(120)

此时取 $c = \sqrt 3 $ ，可使得 $ {w}_{20} = 0，$ 而 ${w}_{0} = 4/9、$ ${w_{10}}$ $= 1/9、$ $ {w}_{11} = 1/36 $ ，即经典的D2Q9模型系数。

对于七阶积分公式，需要满足 $ {H}^{(0)} $ 、 ${H}_{11}^{(2)}$ 、 ${H}_{{1}^{2}{2}^{2}}^{(4)}$ 、 ${H}_{{1}^{4}}^{(4)}$ 、 ${H}_{{1}^{4}{2}^{2}}^{(6)}$ 、 ${H}_{{1}^{6}}^{(6)}$ 对应的6个方程，此时我们尝试用 $ \left|{\xi }_{ai}\right|\leqslant2c $ 积分点组对应的权重，即 ${w_0}～{w_{22}}$ 这6个未知数来求解，其余的未知数设置为0。然而 ${H}_{{1}^{6}}^{(6)}$ 对应的方程为 $ 0 = {r}_{6}，$ 由于 ${r_6}$ 是恒正的，该约束条件在 ${c^2} > 0$ 时无解，故 $ \left|{\xi }_{ai}\right|\leqslant2c $ 无法得到二维情况下的七阶积分公式。事实上，从一维情况下七阶积分必须使用到 ${\xi _a} = \pm 3c$ 的积分点可以推断出这一结论。此时将积分点组扩大到 $ \left|{\xi }_{ai}\right|\leqslant3c $ ，但总未知数个数增加到了10，求解之前需要将部分未知数置零使方程具有定解。由于系数矩阵A的前六行系数组成子矩阵 ${A_{6 \times 10}}$ 的秩为6，我们的目标是寻找 ${A_{6 \times 10}}$ 的六阶非零子式。从总积分点个数最少的原则出发，保留 ${A_{6 \times 10}}$ 中 ${w_0}$ 的列系数（因为它只对应1个积分点），并将 $ {w}_{21}、{w}_{31}、{w}_{32} $ 这三个8积分点组的权重置零，并从 ${A_{6 \times 10}}$ 中移除对应的3列系数；还需从剩下的 ${\boldsymbol{w}_6} = \{$ ${w_{10}},{w_{11}},{w_{20}},{w_{22}},$ ${w_{30}},{w_{33}}\} $ 这6个未知数中选择5个，共有6种可能性。例如，我们将具有最大模 $|{\boldsymbol{\xi }_a}|$ 的权重 ${w_{33}}$ 置零，并移除其列系数，剩下的5个未知数与 ${w_0}$ 对应的列系数构成六阶非零子式，求解该系数子矩阵构成的方程组得：

$ \begin{split}& {w_0} = \frac{{144{c^6} - 392{c^4} + 525{c^2} - 255}}{{144{c^6}}} \\& {w_{10}} = \frac{{36{c^4} - 71{c^2} + 39}}{{48{c^6}}} \\& {w_{11}} = \frac{{4{c^2} - 3}}{{12{c^6}}} \\& {w_{20}} = - \frac{{36{c^4} - 125{c^2} + 75}}{{480{c^6}}} \\& {w_{22}} = \frac{{3 - {c^2}}}{{192{c^6}}} \\& {w_{30}} = \frac{{4{c^4} - 15{c^2} + 15}}{{720{c^6}}} \end{split} $

(121)

此时再选择合适的c使得其中某一个权重为零。由于 ${w_{10}}$ 、 ${w_{30}}$ 都是恒正的， ${w_{11}} = 0$ 对应的 ${c^2} = 3/4$ 与 ${w_{22}} = 0$ 对应的 ${c^2} = 3$ 均会使得某些系数为负值；检验发现只有 ${w_{20}} = 0$ 的一个解可满足所有系数非负的要求，即 ${c^2} = (125 + 5\sqrt {193} )/72$ 。将c代回式（121）即可得到具体的积分点和权重，即Pillippi^[51]中给出的D2V17模型，如图8所示。当然，也可以将 ${{{\boldsymbol{w}}}_6}$ 中 ${w_{33}}$ 以外的一个未知数置零（移除），并移除 ${A_{6 \times 10}}$ 中对应的列向量，执行与上述过程相同的“判定六阶子式非零性→求解方程组→选择合适的c→将c回代得出具体的积分点和权重”这一过程，可求得其他具有17个积分点的七阶积分公式。特别地，将 ${\boldsymbol{w}_6}$ 中 ${w_{20}}$ 或者 ${w_{22}}$ 之一置零（移除），可求解出 ${c^2} = (5 \pm \sqrt {10} )/3$ ，与一维七阶积分公式对应。

对于九阶积分公式，此时需要使用到 ${\boldsymbol{H}^{(0)}}$ ～ ${\boldsymbol{H}^{(8)}}$ 对应的全部8个方程，由于 $H_{{1^8}}^{(8)}$ 对应的方程为 $0 = {r_8}，$ 其解为 ${c^2} \approx 1.19698$ ，与一维九阶积分公式一致。此时未知数 ${w_0}～{w_{33}}$ 有10个，有效的方程只有8个，在保留 ${w_0}$ 的前提下，需要从剩下的9个未知数中选择7个，共有36种选法，其对应的系数矩阵构成八阶子式，再执行与求解七阶积分公式类似的过程。值得注意的是，按最少积分点数原则，若尝试将 ${w_{21}}$ 、 ${w_{31}}$ 、 ${w_{32}}$ 这三个8积分组的未知数取2个置零（移除），由于c的值已经确定，解出的权重总存在一个或者多个负值，违背了权重非负原则。因此，最少积分点序列要保留表2中列出的至少2个8积分点组，其总积分点数最少为 $1 + 4 \times 5 + 8 \times 2 = 37$ 。若保留表2中的前8组积分点组与权重，则得出与Pillippi^[51]相同的D2V37模型，如图9所示。

对于三维情况下的积分公式，也可以遵循同样的思路求解。相关的细节过程可参考文献[47]。

4.4 等温弱可压模型的推导

对于截断到二阶Hermite多项式级数的平衡态分布函数 （式（99）），如果令 $\theta = 1$ ，则有：

$ f_2^{(0)}({{\boldsymbol{\xi}} }) = \rho \omega ({{\boldsymbol{\xi}} })\Big\{ 1 + {\boldsymbol\xi } \cdot {{\boldsymbol{u}}} + \frac{1}{2}\big[{(\boldsymbol{\xi } \cdot {{\boldsymbol{u}}})^2} - {u^2}\big]\Big\} $

(122)

与等温弱可压模型的平衡态分布函数展开式（20）一致。

对于离散速度模型，在4.3节的一维和二维情况已经证实，D1Q3、D2Q9模型的 ${\boldsymbol{\xi }_a}$ 和 ${w_a}$ 分别为一维、二维五阶精度Hermite积分公式的积分点和权重。三维情况下的D3Q15、D3Q19和D3Q27模型也可证明均为三维五阶精度Hermite积分公式的积分点和权重^[47]。

同时，Hermite多项式模型也同时指出（见表1），如需要准确还原动量方程，模拟等温弱可压流动， ${f^{(0)}}$ 至少需要截断到三阶Hermite展开级数，略去 $ {f}_{3}^{（0）} $ 中的三阶量的前提是 $\boldsymbol{u} \to 0$ ，即低马赫数假设，这也与2.1节的结论一致。

5 与传统CFD方法的结合

经典的LBM中时间和空间离散均为沿特征线积分（参考2.2节和4.2节），并设置归一化的 $\Delta t = 1$ ，这极大简化了LBM的计算程序。然而这种时空离散方法也带来了几个弊端：

1）由于 $\Delta x$ 为固定值，只能采用均匀正方形（体）网格。对于边界层流动计算而言，为了解析边界层内的流动而加密全计算域网格会极大增加计算量。在计算量和解析度之间取平衡的办法是使用局部加密网格，但是粗网格和细网格之间的插值会带来新误差^[52]。

2）无法使用贴体网格。在倾斜方向边界，或者曲面边界上，需要同时使用笛卡尔网格和浸没边界法^[53]，以锯齿形边界近似斜边界或者曲面边界。这种非贴体网格由于精度较低，且在对流扩散问题中易造成收敛性困难^[54]和稳定性问题，往往需要与局部加密网格和多松弛时间模型（MRT）结合^[55]。

3）计算稳定性。经验表明，经典LBM方法在雷诺数较大时只能使用较小的物理时间步长才能获得稳定的数值解；当使用Zhou-He非平衡态反弹边界条件^[56]时稳定性问题更加突出，这对实际工程湍流的计算会形成较大的障碍。

由于传统CFD方法可以使用边界层加密贴体网格，使用合适的时间与空间格式时，数值稳定性良好，将传统CFD方法与LBM结合成为了解决以上问题的一种思路。Reider和Sterling^[57]、Cao等^[58]将有限差分法与LBM模型结合，提出了FDLBM方法。FDLBM的基本思路是将LBM方程（式（104））视为关于各个离散分布函数 ${f_a}$ 的双曲型方程，通过有限差分法的时间离散格式，例如Euler格式或者Runger-Kutta格式，将 $\partial {f_a}/\partial t$ 离散化；并使用迎风型或中心型格式将对流项 ${\boldsymbol{\xi }_a} \cdot \nabla {f_a}$ 进行空间离散^[59]。在可压缩问题中，需要将多层速度模型和激波捕捉差分格式结合。Kataoka和Tsutahara^[26]在多层速度模型的基础上，使用二阶迎风格式离散对流项；Wang等^[60]将二阶TVD和五阶WENO格式^[10]与Kataoka-Tsutahara多层速度模型结合使用；许爱国等^[39]则使用了NND格式^[61]来处理对流项，并在一维和二维黎曼问题的计算中取得较为满意的结果。

使用隐式时间离散是传统CFD方法中提高数值稳定性的常用方法，伴随而来的是需要迭代求解线型代数方程组，这往往使得计算量和实施难度显著增加。在FDLBM中，Guo和Zhao^[62]使用了与经典LBM类似的隐式时间处理办法，即类似于式（30）定义一个中间时刻的 ${\bar f_a}$ 代入FDLBM方程中，可将隐式方法转换为显式方法。Wang等^[60]沿用类似的思路，将隐式-显式Runger-Kutta格式^[63]引入FDLBM并消除了隐式格式需要的迭代求解过程。

为了适应复杂几何构型和不规则形状网格的计算需求，Nannell和Succi^[64]、Benzi等^[65]将有限体积方法引入LBM，提出了FVLBM方法。FVLBM的基本思想是，将LBM方程（式（104））对流场中任意一个控制体单元 $\varOmega$ （边界为 $\partial \varOmega$ ，此处 $\varOmega$ 不代表碰撞模型）进行体积分得到：

$ \frac{\partial }{{\partial t}}\int_\varOmega {{f_a}} {\text{d}}\varOmega + \oint_{\partial \varOmega } {({\boldsymbol{\xi }_{a}} \cdot {{\boldsymbol{n}}}){f_a}} {\text{d}}S = - \int_\varOmega {\frac{1}{\tau }} ({f_a} - f_a^{(0)}){\text{d}}\varOmega $

(123)

其中， ${{\boldsymbol{n}}}$ 和 $ {\text{d}}S $ 分别为边界微元 $\partial \varOmega$ 的单位法向量和面积；并定义离散分布函数的单元体平均量：

$ {f_{a,\varOmega }} = \frac{1}{\varOmega }\int_\varOmega {{f_a}} {\text{d}}\varOmega $

(124)

以及单元面对流通量：

$ {L_{a,\varOmega }} = \oint_{\partial \varOmega } {({\boldsymbol{\xi }_{a}} \cdot {{\boldsymbol{n}}}){f_a}} {\text{d}}S $

(125)

则式（123）可写为关于分布函数单元体平均量和单元面对流通量的方程。求解该方程的过程即与传统有限体积方法类似。Chen^[66]和Peng等^[67]将FVLBM的理论和其在非均匀网格上的应用进行了完善。为了提高稳定性，Stiebler等^[68]使用迎风型格式重构单元面通量。Patil和Lakshmisha^[69]则将TVD格式应用于单元面通量的重构。

随着计算流体力学有限元方法日趋成熟^[70-71]，Lee和Lin^[72]、Shi等^[73]将有限元方法引入LBM，发展了有限元LBM方法。近年来，MacMeccan等^[74]、Krüger等^[75]使用混合型有限元LBM方法模拟了可变形颗粒流动并应用于生物流计算领域。由于主题和篇幅所限，此处不对有限元LBM方法作进一步展开。有兴趣的读者可参阅有关文献了解相关细节。

6 总结与展望

LBM多层速度模型起源于单层速度模型，在热流和可压缩流动应用需求下不断发展，经历了平衡态分布函数展开从二阶提高到四阶、离散速度模型从单层网格提高到2层或者更多层网格、约束条件也逐步增多的过程。Hermite多项式模型则提供了一个既理论又实用的判定框架：要想准确还原能量方程，平衡态分布函数必须展开到速度四阶量。Watari-Tsutahara预设的多层速度模型^[24-25]与该结论殊途同归。Hermite多项式模型也同时指出，离散速度模型事实上为相应阶数积分公式的积分点和权重，只需查表（或用程序包计算）找出这些积分点和权重，选出一组应用即可，而不是根据经验预设和优化。尽管Hermite多项式模型给出了较为简洁和确定的形式，但是依照经验理论进行马赫展开或待定系数得到的其他多层速度模型，在特定范围和领域内仍具有一定的应用价值。

必须承认，相对于计算等温弱可压流动的经典单层速度模型而言，多层速度模型仍然是一种尚不成熟的模型，在很多方面仍然需要改进和完善。这主要在以下几个方面：

1）模型的误差、色散与耗散、稳定性分析。目前对单层速度模型的误差和稳定性分析已有不少工作，例如：Sterling和Chen^[76]对D2Q7（六边形模型）、D2Q9和D3Q15进行线性稳定性分析后得出 $\tau > 0.5$ 才能稳定的结论。Marie等^[77]对D3Q19模型进行色散和耗散分析，并与低色散高阶差分格式进行了比较。Silva等^[78]和Bauer等^[79]对最常用的D3Q15、D3Q19和D3Q27模型进行了截断误差分析，指出这三种模型都能在低马赫数下还原动量方程，但是动量输运项的高阶截断误差不同，这将导致在各向异性流动的计算中出现较大差别；在将平衡态分布函数展开式根据连续平衡态分布函数的二阶矩约束条件和D3Q19速度模型进行系数修正后，能获得与D3Q27模型相当的各向同性水平。

对多层速度模型而言，McNamara等^[80]对D3Q21模型（D3Q15加上 $( \pm 2c,0,0)$ 、 $(0, \pm 2c,0)$ 和 $(0,0, \pm 2c)$ 这6个离散点）进行了时间离散格式分析，指出使用Lax-Wendroff时间格式将在高瑞利数下具有更好的稳定性。Wissocq等^[81]则对D2Q9和D2V17进行了模态和稳定性分析，指出D2V17能准确解析线性等温Navier-Stokes方程的3个特征值模态，而D2Q9则在这些模态的高波数上产生一定的误差。这一点与Hermite多项式模型的结论是一致的，即只有使用七阶积分公式的积分点（即D2V17模型）作为离散速度点才能准确还原动量方程；否则由于速度三阶量误差的存在而必须使用低马赫数假设。

值得注意的是，在Hermite多项式模型框架下，经典D3Q15、D3Q19和D3Q27模型从积分精度角度来讲，并不存在差别，因为都具有五阶积分精度；从计算经济性原则来讲，应该选择积分点最少的D3Q15模型。然而Silva等^[78]和Bauer等^[79]的工作提示我们，在选择离散模型时还需要考虑各向同性、截断误差等更多方面的因素，即计算最经济的模型不一定是最优的计算结果。Sieber等^[82]将Hermite多项式模型与D2Q9模型、Chen等^[23]的多层速度模型进行了线性稳定性分析对比，发现对等温弱可压流动而言，将平衡态分布函数展开到三阶Hermite多项式级数，或使用D2V17、D2V37模型都能显著提高稳定性；对于可压缩热流而言，使用四阶Hermite展开的 ${f^{(0)}}$ 与D2V37模型，稳定性都优于Chen等^[23]提出的多层速度模型。该结论在Hermite多项式模型的理论框架下是易于被推断的。然而，对于完全还原了能量方程、且离散速度点数量最少的D2V37、D3V103模型^[47]而言，与同阶积分精度、但具有更多离散速度点的模型横向对比，例如D3V107模型，是否也存在各项同性、稳定性、色散和耗散差异？这一方面需要在后续工作中展开详尽地理论分析；另一方面也需要大量的数值试验，例如方管Poiseuille流动、旋转管道流、强可压缩流动、大温差热流等实际流动计算，来验证这些多层速度模型之间的性质差异。

2）边界条件。LBM单层速度模型在经过多年发展后，其边界条件处理方法已经日趋完善，有关工作可谓汗牛充栋^[33,36,83]，对于壁面反弹边界条件^[56,84]、周期边界条件^[85]、适用于流动出入口的无反射边界条件^[86-87]以及其他类型的边界条件等都已有充分的研究工作，可满足绝大多数计算工况的需求。然而，当使用多层速度模型时，一方面由于离散速度点跨越多层网格，使得边界附近的内部结点与边界结点区分变得模糊；另一方面平衡态分布函数的形式也变得更为复杂，如何恰当处理边界条件成为了一个难点，相关的工作也屈指可数^[83]。Malaspinas等提出^[88]，将边界节点的分布函数离散点按内部和边界分为已知量和未知量，将离散分布函数各阶矩的约束条件构成方程组，通过求解方程组构造一种“通用”的边界条件处理方法。他们在D2Q9和D3Q19模型中进行了测试，并指出该方法可适用于多层速度模型的边界条件（但没有后续验证）。Meng和Zhang^[89]提出了适用于多层速度模型的扩散-反弹边界条件格式，并应用到D2V16、D2V17、D2V37和D3V121^[45]这四种模型，对等温流动和热流进行了测试计算。Lee等^[90]将Hecht和Harting针对D3Q19的非平衡态壁面反弹边界条件^[91]，推广到了用多层速度模型计算等温弱可压流动的情况，并用D2V17模型进行了验算。Klass等^[92]进一步延续了Lee等的工作，并将其应用到了弱可压热流中，使用D2V17和D2V37模型进行了验证，计算结果比Meng和Zhang的扩散-反弹边界条件表现更优。总体来看，针对多层速度模型的边界条件研究工作还处于起步阶段，且大部分局限于弱可压热流和Dirichlet边界类型；对于压缩性较强的跨声速、超声速流动，以及流动出入口涉及的无反射边界条件，尚无相关工作涉足。将单层速度模型中的各类边界条件处理方法扩展到多层速度模型，并在各类流动中广泛验证测试，仍是今后需要较多投入的方向。

3）高性能算法。LBM最大优势是计算速度和高度并行特性，然而它也有一个众所周知的缺点：由于需要存储数量众多的离散分布函数，对计算机内存的使用量比传统CFD方法更多，因此有不少工作聚焦于如何改进算法，减少内存消耗并提高计算效率^[93-95]。另外，正如上一节所述，经典的LBM（相对于FDLBM和FVLBM及其他混合方法而言）一般使用均匀网格，当涉及到边界层计算时，将LBM与自适应局部加密网格（AMR）结合^[96-98]是更加经济有效的方案。这些算法与目前流行的GPU高性能计算相结合，能发挥LBM的最大计算潜力^[99-101]。然而，目前这些高性能算法都是集中应用在单层速度模型中。对于多层速度模型而言，特别是三维情况，准确模拟可压缩热流的最经济模型是D3V103，但其分布函数离散点的数量已经是D3Q27的4倍。这一方面使其内存占用问题更加突出；另一方面，在并行化和局部自适应加密网格时，不同块之间的界面数据交互量也更大，需要相当谨慎地处理。然而这些往往是工程实际应用中决定算法效率的制约因素。发展适用于多层速度模型的高性能算法，包括内存节约算法、高效率通信和AMR方法、CPU+GPU异构高效算法，并在各类复杂几何外形、复杂流动中广泛测试，应引起业内学者的广泛关注。

总体来说，LBM模型，特别是多层速度模型，仍然需要诸多补充和完善才能成为一种可靠、高效的流体力学研究手段和应用技术。正如1998年Chen和Doolen^[2]在文章结尾中提到，LBM多相流和复杂反应系统模型主要集中于等温弱可压流问题，适应于可压缩热流的LBM模型仍待开发。20余年后的今天，这些方面虽然取得了一些进步，但是离可靠性要求还有不少差距。我们期望多层速度模型能够百花齐放，特别是具有较完备数学理论框架的Hermite多项式模型，能在这些方面继续发展壮大并做出应有的贡献。