0 引　言

在过去三十年，格子玻尔兹曼方法逐渐成为计算流体力学领域的一个重要分支，并广泛应用到科学和工程问题当中^[1-2]。不同于基于Navier-Stokes-Fourier（NSF）方程离散的各种传统计算流体力学方法，格子玻尔兹曼方法是一种用统计概率粒子分布函数的碰撞和迁移来描述流体系统的介观数值方法。它与NSF方程的联系可以在小Knudsen数（Kn）下通过Chapman-Enskog（CE）展开^[3]得到。传统计算流体力学方法，如有限差分、有限体积法等，一个难点是处理偏微分方程中的非线性对流项，为了控制数值耗散和数值弥散误差，常采用高阶离散格式。在格子玻尔兹曼方法中，对于通常使用的On-Lattice格子，统计粒子沿特征线进行迁移，即粒子能在一个时间步迁移到达相邻的直角网格格点上。由于没有任何插值操作，仅为计算机上的拷贝操作，这一拉格朗日属性决定了格子玻尔兹曼方法处理对流的简洁性以及它的低数值耗散和数值弥散特性。格子玻尔兹曼方法无需迭代求解压力泊松方程，时间步显式推进，非线性的碰撞项为空间完全局部操作，代数计算为主，数据访存规则，这些特性与目前在高性能计算领域占越来越重要地位的异构并行架构—图形计算处理器（graphic process unit，GPU）架构和国产神威架构（众核处理器）特点相符，为超大规模高精度流体数值模拟提供了可能性^[4-5]。由于格子玻尔兹曼方法的动理学基础、执行简单且高效、天然的并行特性，目前它已被广泛和成功地应用于多相流^[6-7]、湍流^[8-9]、燃烧^[10]、流固耦合^[11-12]等领域。

然而，格子玻尔兹曼方法很多棘手的问题来自于非线性的碰撞项，如数值不稳定性、输运系数不满足伽利略不变性和不可调的普朗特数等。过去几十年，格子玻尔兹曼方法领域的理论工作很多都集中于处理复杂的非线性微分积分碰撞项。最简洁的碰撞模型为钱跃竑等^[13]基于Bhatnagar、Gross和Krook等^[14]提出的Boltzmann-BGK（作者三人首字母命名）方程发展的Lattice BGK(LBGK)模型，该模型中非平衡态部分的各阶成分均以相同的弛豫时间朝平衡态松弛，该算法因执行简单、易于编程而广受欢迎。但是，原始的LBGK模型仅能处理等温流动，数值稳定性存在问题，压力场存在很多虚假噪音。随后很多其他的碰撞模型逐渐被提出，用以实现普朗特数的可调、良好的数值稳定性、虚假压力噪音的压制和消除，以及输运系数的伽利略不变性。对于存在大梯度的流场，几乎所有的计算流体力学数值方法都会遇到数值稳定性问题。对于格子玻尔兹曼方法，在高雷诺数和高马赫数流动中，其数值稳定性问题尤为突出。D'Humières等^[15-16]提出了多松弛时间（multiple-relaxation-time，MRT）碰撞模型。在该方法中，利用正交化手段，通过离散速度集构建了一组正交的基向量，数量与离散速度个数相等。通过这组正交的基向量，粒子分布函数可以计算出不同阶的矩，对不同的矩分配不同的松弛时间，朝对应矩的平衡态松弛。该模型虽然没能通过多松弛因子实现普朗特数的可调，但获得了数值稳定性的大幅提高。另一个能显著提高数值稳定性并改善剪切黏性伽利略不变性的碰撞模型为Geier 等^[17-18]提出的级联多松弛碰撞模型（Cascaded MRT）。所谓的级联多松弛模型即为高阶矩（大于等于三阶）由低阶矩和宏观速度构建，且每阶矩独立松弛。Lycett-Brown等^[19-21]则从一维出发，先通过分布函数与零到二阶矩的关系将分布函数由这些矩显式写出，并用张量积的形式构造了二维和三维级联多松弛碰撞模型。他们的推导相比Geier等的推导更加直观易懂。所谓的级联多松弛模型，本质就是在中心矩空间里面对分布函数的中心矩进行松弛，然后通过中心矩和原点矩的数学关系将中心矩松弛转换到原点矩形式，高阶矩松弛由低阶矩和当地宏观速度表示，即所谓的级联。其他开展级联多松弛模型工作的还有Dubois等^[22-23]、Fei和Luo等^[24-25]，以及De Rosis 等^[26]。然而，由于其采用的离散速度的限制，输运系数的伽利略不变性在级联多松弛模型中并没有完全解决。2015年，Geier 等^[27]受Seeger等^[28]的计算气体动理论累计量方法（Cumulant Method）的启发，提出了累计量格子玻尔兹曼方法（Cumulant LBM），并在后续工作中^[29-30]进行了参数优化讨论。累计量多松弛模型中，以所谓的累计量为独立松弛对象，二阶和三阶矩与前述的级联碰撞模型一样，只有在四阶及以上的矩才产生偏差。Seeger等^[28]在最早的工作中曾指出他的计算气体动理论累积量方法与Grad^[31]等的矩方法在一维空间至少到五阶矩是吻合的，从第六阶矩开始存在偏差。累计量碰撞模型确实取得了良好的数值稳定性，并被应用到湍流边界层模拟^[32]、气动声学模拟^[33]，以及自由面两相流模拟^[34]。Karlin 等^[35-38]提出的熵格子玻尔兹曼模型是另外一个具有良好数值稳定性的碰撞模型。在该模型中，H定理要求在离散分子速度空间得到满足，平衡态分布函数通过迭代求解极值问题得到，高阶格子也通过一维格子的张量积形式得到；松弛因子不再是固定值，而是随流场空间动态变化，类似于大涡模拟模型。熵格子玻尔兹曼模型虽然能显著提高数值稳定性，然而由于其极值问题的迭代求解需要在每个网格点每个时间步执行，其计算效率并无优势。

正则化松弛模型是另外一种具有良好数值稳定性的碰撞模型。最早的正则化碰撞模型由Ladd^[39]于1994年提出，并成功应用到粒子沉降流动中^[40-41]。Ladd认为对于弱可压缩Navier-Stokes（N-S）方程的求解，格子玻尔兹曼碰撞项中只需保留二阶应力张量的松弛（包括剪切应力的松弛和体积应力的松弛），而其他高阶非水动力学动理学模态的松弛频率直接置为1，这样使得不影响N-S方程的非水动力学高阶模态快速松弛。Ladd模型中二阶应力张量中的有迹部分实际为速度空间积分精度误差带来的数值误差项，该项的松弛会调节数值稳定性。实际的二阶应力张量中有迹部分通常只在稠密气体和考虑分子内部自由度激发的时候才会出现，体积黏性松弛频率应置为零，对此Ladd在论文^[41]中做了严格的论述。2005年，Latt等^[42]也提出了一个正则化碰撞模型（类似于Ladd的碰撞模型），并发现正则化执行能大幅提高数值模拟的稳定性和精度；但是他们并没有讨论二阶应力张量中无迹部分和有迹部分的分别松弛，可以视为Ladd模型的一个特殊版本。几乎同时，Chen等^[43]也讨论了与Latt等相同的正则化模型在提高高Kn数流动中数值格式的旋转不变性的前景。Zhang等^[44]利用Shan等^[45]基于Hermite多项式展开Boltzmann-BGK方程的理论体系，将正则化模型延伸到三阶，并在该模型中引入了力项的处理，计算了高Kn数流动。他们认为非平衡态分布函数应该类似于平衡态分布函数，在Hermite截断阶数空间中进行重构；然而该模型所使用的是Hermite原点矩空间，并且在力项的处理中并没有考虑一阶项。2014年，Chen等^[46]将Shan 和Chen在2007年提出的Hermite展开多松弛模型^[47]做了一些改进，将原来的原点矩多松弛改为中心矩多松弛。他们发现二阶矩松弛保持不变，而三阶中心矩松弛转换到原点矩空间后多出了一个修正项，该修正项与局部速度和低阶矩松弛相关，三阶中心矩松弛完全克服了之前模型热传导系数在普朗特数不为1时不满足伽利略不变性的缺陷。2015年，Malaspinas^[48]基于Hermite展开提出了一个递归正则化碰撞模型，证明了非平衡态Hermite系数在等温层级上存在递归关系。该模型中非平衡态Hermite系数的递归关系主要由以下三个关系推导而来：1）平衡态Hermite系数的递归关系，该关系用数学归纳法很容易得到；2）零阶Chapman-Enskog展开得到的欧拉方程；3）一阶Chapman-Enskog展开得到的非平衡态Hermite系数与平衡态Hermite系数的时空导数关系式。该递归正则化碰撞模型被作者证明比原始MRT^[15-16]有更好的数值耗散和色散性质。2017年，Coreixas等^[49]将Malaspinas的模型推广到可压缩传热层级上去。Mattila等^[50]利用中心矩空间Hermite基对非平衡态进行了展开，并在中心矩空间对展开进行了截断，如：对可压缩传热层级的非平衡态展开，将四阶及以上的展开进行截断，相当于将四阶及以上的非平衡态Hermite中心矩设定松弛频率为1；而对等温层级的非平衡态展开，则将三阶及以上的展开进行截断，转换到原点矩空间后发现各阶非平衡态Hermite原点矩刚好与Manaspinas^[48]的递归正则化模型中的非平衡态Hermite矩相同。然而，以上三个正则化模型只对非平衡态Hermite矩进行了讨论，而且只是对各阶非平衡态Hermite原点矩进行独立松弛，并没有在中心矩空间进行松弛，热传导系数的伽利略不变不被保证。2018年，Jacob等^[51]又在Malaspinas的递归正则化模型的基础上提出了一个混合递归正则化模型，对二阶应力张量的重构做了调整；原来正则化模型中该项由非平衡态分布函数的二阶矩得到，在Jacob等的混合正则化模型中，二阶应力张量也可以用速度场的中心差分得到；他们还对两种方法得到的二阶应力张量进行加权，通过调节权值控制额外引入的数值耗散的大小。后续该混合正则化模型被Feng等^[52-54]用到了他们的求解器中。在他们的方法中，质量和动量方程由格子玻尔兹曼方法在D3Q19和D3Q27格子上进行求解，能量方程则通过有限体积方法进行求解，并通过修改平衡态分布函数和添加源项的方式尽可能地消除标准格子存在的速度离散误差，实现了跨声速和超声速模拟。2019年，Shan^[55]在中心矩空间对非平衡态和碰撞项分别进行了展开，并在中心矩空间对每阶矩进行独立松弛，然后通过中心矩空间和原点矩空间Hermite基之间的关系将碰撞项转换回原点矩空间，得到了递归形式的碰撞项。该碰撞模型继承了2007年Shan等^[47]提出的原点矩空间Hermite矩多松弛模型的优点，克服了Chen等^[46]提到的热传导系数的伽利略不变性问题，并可以推广到任意高阶松弛。随后，Li等^[56]在温度标度的中心矩空间中，对非平衡态和碰撞项进行了展开，并在该空间中对各阶矩进行松弛，利用温度标度中心矩空间和原点矩空间Hermite基之间的关系将碰撞项转换回原点矩空间，也得到了递归形式的碰撞项；与原点矩空间松弛的碰撞项^[47]相比，三阶及以上的碰撞项包含速度、温度及低阶矩松弛构成的修正项。

关于正则化碰撞模型的边界条件处理，历史上也有一些学者取得了部分进展。2007年，Niu等^[57]模拟了微尺度稀薄气体流动，讨论了散射边界条件、松弛时间和正则化碰撞模型。2008年，Latt等^[58]基于正则化重构非平衡态的思路发展了适用于平直边界的正则化边界条件。2011年，Malaspinas等^[59]提出了一个通用的处理平直边界的正则化边界条件，通过建立边界点已知分布函数与宏观量的超定方程组，解出宏观矩，再重构出边界点的非平衡态分布函数。该边界处理方法可适用于标准格子（D2Q9、D3Q19等），也适用于高阶格子（D2V37、D3V103等），但是否适用于任意曲线曲面边界需进一步研究。2017年，Wissocq等^[60]发展了适用于声学模拟的正则化特征边界条件。2019年，Guo等^[61]在混合正则化模型中讨论了可压缩流动的任意曲面边界条件和域边界特征边界条件。

在将正则化碰撞模型应用到多相流方面，也有一些进展。Otomo等将正则化碰撞模型分别应用到多组分伪势模型^[62]和相场多相流模型^[63]。Ba等^[64]将正则化碰撞模型应用到颜色梯度模型。另外，将正则化碰撞模型用于噪声模拟也有很多代表性工作，如Zhuo等^[65]、Chen等^[66]、Casalino等^[67-68]、Gendre等^[69]、Astoul等^[70]。本文对正则化碰撞模型在各个领域的应用不再做详细展开。

本文将重点阐述、回顾正则化碰撞模型，通过一个系统的理论分析框架阐述作者与合作者提出的高阶正则化碰撞模型，并对以往的正则化碰撞模型的理论关系进行梳理，同时对正则化多松弛模型的本质进行归纳分析。后文将按以下顺序进行论述：第1节阐述Hermite展开的基础理论框架，为后文中的碰撞算子正则化提供分析的准则；第2节则利用第1节建立的理论体系，对正则化碰撞模型发展历史上有重要贡献的经典模型逐个进行分析、梳理和比较；最后对全文进行归纳总结，得出全文理论分析的结论。

1 理论基础

该部分我们从Boltzmann方程出发，在Hermite谱空间对分布函数进行展开，分别在原点矩空间、中心矩空间和温度标度的中心矩空间对Maxwell-Boltzmann平衡态、非平衡态、碰撞项和力项进行展开，为后面的正则化理论模型回顾做准备。

我们先从Boltzmann方程进行阐述。Boltzmann方程为基于分子二体碰撞、分子混沌假设和外力不影响局部碰撞的动力学行为等三个假设，推导的稀薄气体系统的控制方程。它描述了统计粒子分布函数在速度空间和物理空间的演化，其右端的非线性碰撞项极为复杂，直接求解困难重重。Boltzmann-BGK方程^[14]作为Boltzmann方程碰撞项的一个极度简化版本，认为所有动理学矩都以相同的速度朝平衡态弛豫，其控制方程如下：

$ \frac{{\partial f}}{{\partial t}} + {\boldsymbol{{{{\boldsymbol{\xi}} }} }} \cdot \nabla f + {\boldsymbol{g}} \cdot {\nabla _{\boldsymbol{{{{\boldsymbol{\xi}} }} }}}f = \varOmega = - \frac{1}{\tau }(f - {f^{{\rm{eq}}}}) $

(1)

其中 $ f({\boldsymbol{x}},{\boldsymbol{{{{\boldsymbol{\xi}} }} }},t) $ 代表统计粒子在特定的分子速度 $ {\boldsymbol{{{{\boldsymbol{\xi}} }} }} $ 、物理空间位置 $ {\boldsymbol{x}} $ 和时间 $ t $ 的概率。式（1） 左边为统计粒子的时空输运，右边为粒子的局部碰撞。选择特征速度为 $ {c_s} = \sqrt {R{T_0}/{m_0}} $ （其中， $ R $ 为Boltzmann 常数， $ {T_0} $ 为参考温度， $ {m_0} $ 为气体分子的单位质量），式（1）中的Maxwell-Boltzmann 平衡态分布函数为：

$ {f^{{\rm{eq}}}}({\boldsymbol{x}},{\boldsymbol{{{{\boldsymbol{\xi}} }} }},t) = \frac{\rho }{{{{(2\text{π} \theta )}^{D/2}}}}\exp \Bigg[ - \frac{{{{({\boldsymbol{{{{\boldsymbol{\xi}} }} }} - {\boldsymbol{u}})}^2}}}{{2\theta }}\Bigg] $

(2)

其中， $ \theta $ 为无量纲温度，D为空间维度。如果 $\; {\rho _0} $ 设定为参考密度，则无量纲粒子分布函数可定义为 $ \bar f = f{c_s}^D/{\rho _0} $ 。为简便起见，在下文中省略分布函数上的一拔标记。

1.1 粒子分布函数的Hermite展开

因为Hermite多项式良好的正交特性，Grad^[31]在1949年早期的工作中已将其用到稀薄气体动理学理论中。本论文中，我们分别使用三种不同自变量的Hermite基 $ {\boldsymbol{H}}({\boldsymbol{{{{\boldsymbol{\xi}} }} }}) $ 、 $ {\boldsymbol{H}}({\boldsymbol{c}}) $ 、 $ {\boldsymbol{H}}({\boldsymbol{\eta }}) $ （其中 $ {\boldsymbol{c}} = {\boldsymbol{{{{\boldsymbol{\xi}} }} }} - {\boldsymbol{u}} $ ， $ {\boldsymbol{\eta }} = {{\boldsymbol{c}} \mathord{\left/ {\vphantom {{\boldsymbol{c}} {\sqrt \theta }}} \right. } {\sqrt \theta }} $ ）来展开分布函数：第一个基的自变量为绝对坐标系中的分子速度；第二个基的自变量为平移坐标系中的分子速度；第三个基的自变量为平移坐标系中的局部马赫数（仅差一个绝热指数），该自变量为无量纲量。粒子分布函数的Hermite展开可视为一种谱展开。原点矩空间中粒子分布函数可展开如下：

$ f({\boldsymbol{x}},{\boldsymbol{{{{\boldsymbol{\xi}} }} }},t) = \omega ({\boldsymbol{{{{\boldsymbol{\xi}} }} }})\sum\limits_{n = 0}^\infty {\frac{1}{{n!}}{{\boldsymbol{a}}^{(n)}}:{{\boldsymbol{H}}^{(n)}}({\boldsymbol{{{{\boldsymbol{\xi}} }} }})} $

(3)

上式中右端项的符号"："代表两个n阶张量的全缩并（下文同）。

原点矩Hermite基定义如下：

$ {{\boldsymbol{H}}^{(n)}}({\boldsymbol{{{{\boldsymbol{\xi}} }} }}) = \frac{{{{( - 1)}^n}}}{{\omega ({\boldsymbol{{{{\boldsymbol{\xi}} }} }})}}{\nabla ^n}\omega ({\boldsymbol{{{{\boldsymbol{\xi}} }} }}) $

(4)

相应的权函数为：

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

(5)

其中 $ {{{{\boldsymbol{\xi}} }} ^2} = {\boldsymbol{{{{\boldsymbol{\xi}} }} }} \cdot {\boldsymbol{{{{\boldsymbol{\xi}} }} }} $ 。

式（3）两端同乘以 $ {{\boldsymbol{H}}^{(n)}}({\boldsymbol{{{{\boldsymbol{\xi}} }} }}) $ ，并在分子速度空间中积分，利用Hermite多项式的正交特性，可以求解相应阶的Hermite系数为：

$ {{\boldsymbol{a}}^{(n)}} = \int {f{\boldsymbol{H}^{(n)}}({\boldsymbol{{{{\boldsymbol{\xi}} }} }}){\rm{d}}{\boldsymbol{{{{\boldsymbol{\xi}} }} }}} $

(6)

1.2 Hermite展开的截断和速度空间离散

在实际的数值执行中，粒子分布函数不可能展开到无穷阶，必须进行截断。截断的阶数直接影响水动力学宏观方程的适用范围及其精度，这意味着实际上需要将分布函数投影到有限阶Hermite基上去。对于截断阶数如何影响恢复的水动力学方程层级以及其精度，Shan等^[45]利用Chapman-Enskog（CE）对其进行了深入分析，这里简略地回顾一下。设 $ \varepsilon $ 为与Kn数同量级的小量，分布函数以及时间和空间导数可以按 $ \varepsilon $ 量级进行摄动展开，如下：

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

(7)

$ {\partial _t} = \varepsilon \partial _t^{(0)} + {\varepsilon ^2}\partial _t^{(1)} + \cdots $

(8)

$ \nabla = \varepsilon \nabla $

(9)

将上述几式代入Boltzman-BGK方程，并将 $ \varepsilon $ 同阶量整理出来，可得到下面k+1阶分布函数与k阶及以下分布函数的时空导数递归关系：

$ {f^{(k + 1)}} = - \tau \left(\sum\limits_{m = 0}^k {\partial _t^{(m)}{f^{(k - m)}} + {\boldsymbol{{{{\boldsymbol{\xi}} }} }}} \cdot \nabla {f^{(k)}}\right) $

(10)

例如，一阶和二阶递归关系可显式写出：

$ {f^{(1)}} = - \tau \Big(\partial _t^{(0)} + {\boldsymbol{{{{\boldsymbol{\xi}} }} }} \cdot \nabla \Big){f^{(0)}} $

(11)

$ {f^{(2)}} = - \tau \bigg[\Big(\partial _t^{(0)} + {\boldsymbol{{{{\boldsymbol{\xi}} }} }} \cdot \nabla \Big){f^{(0)}} + \partial _t^{(1)}{f^{(1)}}\bigg] $

(12)

类似地，k阶分布函数亦可用Hermite级数展开如下：

$ {f^{(k)}} = \omega ({\boldsymbol{{{{\boldsymbol{\xi}} }} }})\sum\limits_{n = 0}^\infty {\frac{1}{{n!}}{\boldsymbol{a}}_k^{(n)}:{{\boldsymbol{H}}^{(n)}}({\boldsymbol{{{{\boldsymbol{\xi}} }} }})} $

(13)

利用递归关系式：

$ {{\boldsymbol{\xi }}_i}{{\boldsymbol{H}}^{(n)}}({\boldsymbol{\xi }}) = {{\boldsymbol{H}}^{(n + 1)}}({\boldsymbol{\xi }}) + {{\bf{\delta }}_i}{{\boldsymbol{H}}^{(n - 1)}}({\boldsymbol{\xi }}) $

(14)

以及Hermite多项式的正交性，可得到下式：

$ \int {{{\boldsymbol{H}}^{(n)}}({\boldsymbol{\xi }})({\boldsymbol{\xi }} \cdot \nabla f)d{\boldsymbol{\xi }}} = \nabla \cdot {{\boldsymbol{a}}^{(n + 1)}} + n\nabla {{\boldsymbol{a}}^{(n - 1)}}$

(15)

将式（13）代入式（10），并利用式（15），则可得到如下Hermite系数之间的递归关系式：

$ {\boldsymbol{a}}_{k + 1}^{(n)} = - \tau \Bigg(\sum\limits_{m = 0}^k {\partial _t^{(m)}{\boldsymbol{a}}_{k - m}^{(n)} + n\nabla {\boldsymbol{a}}_k^{(n - 1)} + } \nabla \cdot {\boldsymbol{a}}_k^{(n + 1)}\Bigg) $

(16)

这里要注意一点，k代表CE展开中的阶数，n代表Hermite展开中的阶数，每一阶的CE展开中，都有完整的Hermite展开。从上式可以看出，k+1阶CE展开中n阶Hermite系数只依赖于k阶及以下CE展开中最多n+1阶Hermite系数。这意味着如果将平衡态分布函数 $ {f^{(0)}} $ 的Hermite展开截断到n阶，则k阶CE展开中前n-k阶Hermite系数与使用平衡态分布函数的全展开表达式是等价的。通常所求解的NSF水动力学宏观方程需要涉及二阶应力张量和三阶热流矢量，对应非平衡态分布函数量 $ {f^{(1)}} $ 的二阶和三阶Hermite系数。这意味着如果要精确表达二阶应力张量和三阶热流矢量，则平衡态的Hermite展开要分别至少截断到三阶和四阶。因此，后续的Hermite展开我们均不超过四阶。Hermite系数的递归关系式（16）在后面正则化模型分析中将会用到。

受益于Hermite多项式的正交属性，粒子分布函数的前若干阶矩是可以由截断形式的Hermite展开得到的粒子分布函数获得。因此，粒子分布函数可以投影到由前N阶Hermite多项式构成的Hilbert空间，并保证其前N阶矩与完整粒子分布函数的前N阶矩严格相等。粒子分布函数的前N阶Hermite展开截断为：

$ f({\boldsymbol{x}},{\boldsymbol{{{{\boldsymbol{\xi}} }} }},t) \approx {f^N}({\boldsymbol{x}},{\boldsymbol{{{{\boldsymbol{\xi}} }} }},t) = \omega ({\boldsymbol{{{{\boldsymbol{\xi}} }} }})\sum\limits_{n = 0}^N {\frac{1}{{n!}}{{\boldsymbol{a}}^{(n)}}({\boldsymbol{x}},t):{{\boldsymbol{H}}^{(n)}}({\boldsymbol{{{{\boldsymbol{\xi}} }} }})} $

(17)

根据Hermite基的正交特性，投影到无穷阶Hermite基上和投影到截断至N阶的Hermite基上，要使两种投影方式下M阶Hermite系数相等，必须满足 $ M \leqslant N $ ，即：

$ {{\boldsymbol{a}}^{(M)}}({\boldsymbol{x}},t) = \int {f({\boldsymbol{x}},{\boldsymbol{{{{\boldsymbol{\xi}} }} }},t){{\boldsymbol{H}}^{(M)}}({\boldsymbol{{{{\boldsymbol{\xi}} }} }}){\rm{d}}{\boldsymbol{{{{\boldsymbol{\xi}} }} }}} = \int {{f^N}({\boldsymbol{x}},{\boldsymbol{{{{\boldsymbol{\xi}} }} }},t){{\boldsymbol{H}}^{(M)}}({\boldsymbol{{{{\boldsymbol{\xi}} }} }}){\rm{d}}{\boldsymbol{{{{\boldsymbol{\xi}} }} }}} $

(18)

上式右端积分的被积函数为M+N的多项式，积分可以在离散速度空间上选定特定的速度集^{[45, 71]}，用Gauss-Hermite积分规则可精确算出，即：

$ \int {\omega ({\boldsymbol{{{{\boldsymbol{\xi}} }} }})\frac{{{f^N}({\boldsymbol{x}},{\boldsymbol{{{{\boldsymbol{\xi}} }} }},t)}}{{\omega ({\boldsymbol{{{{\boldsymbol{\xi}} }} }})}}{{\boldsymbol{H}}^{(M)}}({\boldsymbol{{{{\boldsymbol{\xi}} }} }}){\rm{d}}{\boldsymbol{{{{\boldsymbol{\xi}} }} }}} = \sum\limits_0^{b - 1} {{w_i}\frac{{{f^N}({\boldsymbol{x}},{{\boldsymbol{{{{\boldsymbol{\xi}} }} }}_i},t)}}{{\omega ({{\boldsymbol{{{{\boldsymbol{\xi}} }} }}_i})}}{{\boldsymbol{H}}^{(M)}}({{\boldsymbol{{{{\boldsymbol{\xi}} }} }}_i})} $

(19)

上式中等号严格成立的前提是Gauss-Hermite积分精度Q满足不等式 $ M + N \leqslant Q $ 。速度相空间离散形式的粒子分布函数可以定义为 $ {f_i} = {w_i}{{{f^N}({\boldsymbol{x}},{{\boldsymbol{{{{\boldsymbol{\xi}} }} }}_i},t)} \mathord{\left/ {\vphantom {{{f^N}({\boldsymbol{x}},{{\boldsymbol{{{{\boldsymbol{\xi}} }} }}_i},t)} {\omega ({{\boldsymbol{{{{\boldsymbol{\xi}} }} }}_i})}}} \right. } {\omega ({{\boldsymbol{{{{\boldsymbol{\xi}} }} }}_i})}} $ 。因此，Hermite展开截断形式和速度空间离散形式的分布函数的矩与连续分布函数的矩等价。数学条件为矩的阶数小于等于Hermite展开的阶数，Gauss-Hermite积分精度大于等于二者之和^{[45, 71]}，即：

$ M \leqslant N \;\;和 \;\; M + N \leqslant Q $

(20)

1.3 不同Hermite基之间的数学联系

为了便于后文对不同空间中Hermite系数进行评估，对于不同Hermite基之间的转换关系做一个简单介绍，详细的推导过程可以参考Shan^{[55, 72]}、Li等^[56]、Shan等^[72]的论文附录部分，他们给出了详细的数学推导，这里只给出前四阶的数学转换关系。由于在实际的数值执行中，中心矩空间的离散速度和权值会随当地的宏观速度而变化；而温度标度中心矩空间的离散速度和权值则会随当地的宏观速度和温度发生变化。这会破坏格子玻尔兹曼方法完美的迁移操作，而引入额外的数值耗散和色散，让格子玻尔兹曼方法的优点丧失。因此，我们的思路是在中心矩空间和温度标度中心矩空间完成碰撞，然后通过它们与原点矩空间的数学关系转换到原点矩空间，最终在原点矩空间实现粒子迁移操作。

首先写出原点矩Hermite基的零到四阶表达式，如下（具体可参考文献[55]）：

$\begin{split}& {{\boldsymbol{H}}^{(0)}}({\boldsymbol{{{{\boldsymbol{\xi}} }} }}) = 1\\& {{\boldsymbol{H}}^{(1)}}({\boldsymbol{{{{\boldsymbol{\xi}} }} }}) = {\boldsymbol{{{{\boldsymbol{\xi}} }} }} \\& {{\boldsymbol{H}}^{(2)}}({\boldsymbol{{{{\boldsymbol{\xi}} }} }}) = {\boldsymbol{{{{\boldsymbol{\xi}} }} {{{\boldsymbol{\xi}} }} }} - {\boldsymbol{\delta }} \\& {{\boldsymbol{H}}^{(3)}}({{{{{\boldsymbol{\xi}} }} }}) = {{{{{\boldsymbol{\xi}} }} {{{\boldsymbol{\xi}} }} {{{\boldsymbol{\xi}} }} }} - 3{{{{{\boldsymbol{\xi}} }}\boldsymbol\delta }}\\& {{\boldsymbol{H}}^{(4)}}({\boldsymbol{{{{\boldsymbol{\xi}} }} }}) = {\boldsymbol{{{{\boldsymbol{\xi}} }} {{{\boldsymbol{\xi}} }} {{{\boldsymbol{\xi}} }} {{{\boldsymbol{\xi}} }} }} - 6{\boldsymbol{{{{\boldsymbol{\xi}} }} {{{\boldsymbol{\xi}} }} \delta }} + 3{\boldsymbol{\delta \delta }}\\\end{split}$

(21)

中心矩Hermite基和原点矩Hermite基之间的数学关系可通过二项式转换得到^{[49, 55]}，写出零到四阶的显式转换关系如下：

$\begin{split}&{{\boldsymbol{H}}^{(0)}}({\boldsymbol{c}}) = 1\\& {{\boldsymbol{H}}^{(1)}}({\boldsymbol{c}}) = {{\boldsymbol{H}}^{(1)}}({\boldsymbol{{{{\boldsymbol{\xi}} }} }}) - {\boldsymbol{u}}\\& {{\boldsymbol{H}}^{(2)}}({\boldsymbol{c}}) = {{\boldsymbol{H}}^{(2)}}({\boldsymbol{{{{\boldsymbol{\xi}} }} }}) - 2{{\boldsymbol{H}}^{(1)}}({\boldsymbol{{{{\boldsymbol{\xi}} }} }}){\boldsymbol{u}} + {\boldsymbol{uu}}\\& {{\boldsymbol{H}}^{(3)}}({\boldsymbol{c}}) = {{\boldsymbol{H}}^{(3)}}({\boldsymbol{{{{\boldsymbol{\xi}} }} }}) - 3{{\boldsymbol{H}}^{(2)}}({\boldsymbol{{{{\boldsymbol{\xi}} }} }}){\boldsymbol{u}} + 3{{\boldsymbol{H}}^{(1)}}({\boldsymbol{{{{\boldsymbol{\xi}} }} }}){\boldsymbol{uu}} - {\boldsymbol{uuu}} \\& {{\boldsymbol{H}}^{(4)}}({\boldsymbol{c}}) = {{\boldsymbol{H}}^{(4)}}({\boldsymbol{{{{\boldsymbol{\xi}} }} }}) - 4{{\boldsymbol{H}}^{(3)}}({\boldsymbol{{{{\boldsymbol{\xi}} }} }}){\boldsymbol{u}} + 6{{\boldsymbol{H}}^{(2)}}({\boldsymbol{{{{\boldsymbol{\xi}} }} }}){\boldsymbol{uu}} -\\& \qquad\qquad 4{{\boldsymbol{H}}^{(1)}}({\boldsymbol{{{{\boldsymbol{\xi}} }} }}){\boldsymbol{uuu}} + {\boldsymbol{uuuu}}\\ \end{split}$

(22)

温度标度中心矩Hermite基和原点矩Hermite基之间的数学关系可通过数学归纳法证明得到^{[56, 72]}，写出零到四阶的显式转换关系如下：

$\begin{split}&{{\boldsymbol{H}}^{(0)}}({\boldsymbol{\eta }}) = 1\\& {{\boldsymbol{H}}^{(1)}}({\boldsymbol{\eta }}) = {\theta ^{ - \frac{1}{2}}}\Big[{{\boldsymbol{H}}^{(1)}}({\boldsymbol{{{{\boldsymbol{\xi}} }} }}) - {\boldsymbol{u}}\Big]\\& {{\boldsymbol{H}}^{(2)}}({\boldsymbol{\eta }}) = {\theta ^{ - \frac{2}{2}}}\Big[{{\boldsymbol{H}}^{(2)}}({\boldsymbol{{{{\boldsymbol{\xi}} }} }}) - 2{{\boldsymbol{H}}^{(1)}}({\boldsymbol{{{{\boldsymbol{\xi}} }} }}){\boldsymbol{u}} + {{\boldsymbol{u}}_1}\Big]\\& {{\boldsymbol{H}}^{(3)}}({\boldsymbol{\eta }}) = {\theta ^{ - \frac{3}{2}}}\Big[{{\boldsymbol{H}}^{(3)}}({\boldsymbol{{{{\boldsymbol{\xi}} }} }}) - 3{{\boldsymbol{H}}^{(2)}}({\boldsymbol{{{{\boldsymbol{\xi}} }} }}){\boldsymbol{u}} + 3{{\boldsymbol{H}}^{(1)}}({\boldsymbol{{{{\boldsymbol{\xi}} }} }}){{\boldsymbol{u}}_1} - {\boldsymbol{u}}{{\boldsymbol{u}}_2}\Big]\\& {{\boldsymbol{H}}^{(4)}}({\boldsymbol{\eta }}) = {\theta ^{ - \frac{4}{2}}}\Big[{{\boldsymbol{H}}^{(4)}}({\boldsymbol{{{{\boldsymbol{\xi}} }} }}) - 4{{\boldsymbol{H}}^{(3)}}({\boldsymbol{{{{\boldsymbol{\xi}} }} }}){\boldsymbol{u}} + 6{{\boldsymbol{H}}^{(2)}}({\boldsymbol{{{{\boldsymbol{\xi}} }} }}){{\boldsymbol{u}}_1}-\\&\qquad\qquad 4{{\boldsymbol{H}}^{(1)}}({\boldsymbol{{{{\boldsymbol{\xi}} }} }}){\boldsymbol{u}}{{\boldsymbol{u}}_2} + {\boldsymbol{uu}}{{\boldsymbol{u}}_2} - 3(\theta - 1){{\boldsymbol{u}}_1}{\boldsymbol{\delta }}\Big]\\ \end{split}$

(23)

其中， ${{\boldsymbol{u}}_1} \equiv {{\boldsymbol{u}}^2} - \left( {\theta - 1} \right){\boldsymbol{\delta}}$ ， ${{\boldsymbol{u}}_2} \equiv {{\boldsymbol{u}}^2} - 3\left( {\theta - 1} \right){\boldsymbol{\delta}}$ 。式（21）～式（23）以及后文中的高阶对称张量，其前面的系数表示轮换的项数，如 $ 4{{\boldsymbol{H}}^{(1)}}({\boldsymbol{{{{\boldsymbol{\xi}} }} }}){\boldsymbol{uuu}} = {H_{\alpha} }{u_\beta }{u_\gamma }{u_\delta } + {H_\beta }{u_\gamma }{u_\delta }{u_{\alpha} } + {H_\gamma }{u_\delta }{u_{\alpha} }{u_\beta } + {H_\delta }{u_{\alpha} }{u_\beta }{u_\gamma } $ ，而其中的下标 ${\alpha} 、\beta 、\gamma 、\delta$ 可取三维空间坐标分量 $x、y、z$ 。

1.4 不同Hermite基级数展开及其系数评估

实际上，历史上发展的不同的正则化模型的差异主要在以下几个地方：一个是使用的Hermite基不同，即在不同的空间（原点矩空间、中心矩空间和温度标度中心矩空间）中展开；另一个是分布函数Hermite展开截断的阶数的差异，尤其是非平衡态Hermite展开的截断；还有一个是碰撞项展开的空间和碰撞项中弛豫的张量矩分量。因此，为了对这些正则化碰撞模型进行系统地回顾，我们将平衡态、非平衡态、碰撞项和力项都用三种不同的Hermite基进行展开，然后对它们进行系统地分析。

1.4.1 平衡态分布函数的Hermite展开

使用权函数定义表达式（5），Maxwell-Boltzmann 平衡态分布函数可写为下式：

$ {f^{(0)}} = \frac{\rho }{{{{\Big(\sqrt \theta \Big)}^D}}}\omega \bigg(\frac{{{\boldsymbol{{{{\boldsymbol{\xi}} }} }} - {\boldsymbol{u}}}}{{\sqrt \theta }}\bigg) $

(24)

类似地，平衡态分布函数可用它在前N阶Hermite多项式构成的Hilbert空间上的投影来近似，即：

$ {f^{(0)}}({\boldsymbol{x}},{\boldsymbol{{{{\boldsymbol{\xi}} }} }},t) \approx {f^{(0),N}}({\boldsymbol{x}},{\boldsymbol{{{{\boldsymbol{\xi}} }} }},t) = \omega ({\boldsymbol{{{{\boldsymbol{\xi}} }} }})\sum\limits_{n = 0}^N {\frac{1}{{n!}}{\boldsymbol{a}}_0^{(n)}:{{\boldsymbol{H}}^{(n)}}({\boldsymbol{{{{\boldsymbol{\xi}} }} }})} $

(25)

则平衡态分布函数的原点矩Hermite系数可由下式计算得到：

$ {\boldsymbol{a}}_0^{(n)} = \int {{f^{(0)}}{{\boldsymbol{H}}^{(n)}}({\boldsymbol{{{{\boldsymbol{\xi}} }} }}){\rm{d}}{\boldsymbol{{{{\boldsymbol{\xi}} }} }}} $

(26)

将Maxwell-Boltzmann平衡态分布函数代入上式即可得到相应阶的原点矩Hermite系数^[45]，显式地写出零到四阶有：

$\begin{split}&{\boldsymbol{a}}_0^{(0)} = \rho\\& {\boldsymbol{a}}_0^{(1)} = \rho {\boldsymbol{u}}\\& {\boldsymbol{a}}_0^{(2)} = \rho \Big[{{\boldsymbol{u}}^2} + (\theta - 1){\boldsymbol{\delta }}\Big]\\& {\boldsymbol{a}}_0^{(3)} = \rho \Big[{{\boldsymbol{u}}^3} + 3(\theta - 1){\boldsymbol{\delta u}}\Big]\\& {\boldsymbol{a}}_0^{(4)} = \rho \Big[{{\boldsymbol{u}}^4} + 6(\theta - 1){\boldsymbol{\delta }}{{\boldsymbol{u}}^2} + 3{(\theta - 1)^2}{{\boldsymbol{\delta }}^2}\Big]\\ \end{split}$

(27)

如果平衡态分布函数在中心矩空间进行Hermite展开，即用Hermite基 $ {{\boldsymbol{H}}^{\left( n \right)}}({\boldsymbol{{{{\boldsymbol{\xi}} }} }} - {\boldsymbol{u}}) $ 进行级数展开：

$ {f^{(0)}}({\boldsymbol{x}},{\boldsymbol{{{{\boldsymbol{\xi}} }} }},t) \approx {f^{(0),N}}({\boldsymbol{x}},{\boldsymbol{{{{\boldsymbol{\xi}} }} }},t) = \omega ({\boldsymbol{c}})\sum\limits_{n = 0}^N {\frac{1}{{n!}}{\boldsymbol{b}}_0^{(n)}:{{\boldsymbol{H}}^{(n)}}({\boldsymbol{c}})} $

(28)

其中 $ {\boldsymbol{c}} = {\boldsymbol{{{{\boldsymbol{\xi}} }} }} - {\boldsymbol{u}} $ 。相应的中心矩Hermite系数亦可根据正交性得到：

$ {\boldsymbol{b}}_0^{(n)} = \int {{f^{(0)}}{{\boldsymbol{H}}^{(n)}}({\boldsymbol{c}})d{\boldsymbol{c}}} $

(29)

利用中心矩Hermite基和原点矩Hermite基的数学转换关系式（22），上式中各阶平衡态中心矩Hermite系数可相应得到，这里显式地写出零到四阶：

$\begin{split}&{\boldsymbol{b}}_0^{(0)} = \rho\\& {\boldsymbol{b}}_0^{(1)} = 0\\& {\boldsymbol{b}}_0^{(2)} = \rho (\theta - 1)\boldsymbol\delta\\& {\boldsymbol{b}}_0^{(3)} = 0\\& {\boldsymbol{b}}_0^{(4)} = 3\rho {(\theta - 1)^2}{\boldsymbol{\delta \delta }}\\\end{split} $

(30)

对于平衡态的中心矩Hermite系数，奇数阶为零，且系数中不再出现局部宏观速度，即各阶矩与坐标系无关。

如果平衡态分布函数在温度标度的中心矩空间进行展开，即基于Hermite基 $ {{\boldsymbol{H}}^{\left( n \right)}}({{\boldsymbol{\eta}} }) $ 对平衡态进行级数展开，则有：

${f^{(0)}}({\boldsymbol{x}},{\boldsymbol{{{{\boldsymbol{\xi}} }} }},t) \approx {f^{(0),N}}({\boldsymbol{x}},{\boldsymbol{{{{\boldsymbol{\xi}} }} }},t) = \omega ({\boldsymbol{\eta }})\sum\limits_{n = 0}^N {\frac{1}{{n!}}{\boldsymbol{d}}_0^{(n)}:{{\boldsymbol{H}}^{(n)}}({{\boldsymbol{\eta}}})} $

(31)

其中 $ {\boldsymbol{\eta }} = {{({\boldsymbol{{{{\boldsymbol{\xi}} }} }} - {\boldsymbol{u}})} \mathord{\left/ {\vphantom {{({\boldsymbol{{{{\boldsymbol{\xi}} }} }} - {\boldsymbol{u}})} {\sqrt \theta }}} \right. } {\sqrt \theta }} $ 。相应地，平衡态分布函数的各阶温度标度中心矩Hermite系数为：

$ {\boldsymbol{d}}_0^{(n)} = \int {{f^{(0)}}{{\boldsymbol{H}}^{(n)}}({\boldsymbol{\eta }})\frac{1}{{{{\Big(\sqrt \theta \Big)}^D}}}{\rm{d}}{\boldsymbol{{{{\boldsymbol{\xi}} }} }}} $

(32)

上式中 ${\rm{d}}{{\boldsymbol{\eta}}} = {{{\rm{d}}{\boldsymbol{{{{\boldsymbol{\xi}} }} }}} \mathord{\left/ {\vphantom {{d{\boldsymbol{{{{\boldsymbol{\xi}} }} }}} {{{(\sqrt \theta )}^D}}}} \right. } {{{\Big(\sqrt \theta \Big)}^D}}}$ 。利用式（24）中平衡态的表达式，以及Hermite多项式的正交特性可知，上式积分只有零阶项不为零，其他项均为零，即：

$ {\boldsymbol{d}}_0^{(0)} = \rho /{\Big(\sqrt \theta \Big)^D} $

(33)

1.4.2 非平衡态分布函数的Hermite展开

类似地，非平衡态分布函数亦可用前N阶Hermite多项式构成的Hilbert空间上的投影来近似。如果选择原点矩Hermite基进行展开，则有：

$ {f^{({{1}})}}({\boldsymbol{x}},{\boldsymbol{{{{\boldsymbol{\xi}} }} }},t) \approx {f^{(1),N}}({\boldsymbol{x}},{\boldsymbol{{{{\boldsymbol{\xi}} }} }},t) = \omega ({\boldsymbol{{{{\boldsymbol{\xi}} }} }})\sum\limits_{n = 0}^N {\frac{1}{{n!}}{\boldsymbol{a}}_1^{(n)}:{{\boldsymbol{H}}^{(n)}}({\boldsymbol{{{{\boldsymbol{\xi}} }} }})} $

(34)

相应地，非平衡态原点矩Hermite系数可由下式计算：

$ {\boldsymbol{a}}_1^{(n)} = \int {{f^{(1)}}{{\boldsymbol{H}}^{(n)}}({\boldsymbol{{{{\boldsymbol{\xi}} }} }}){\rm{d}}{\boldsymbol{{{{\boldsymbol{\xi}} }} }}} $

(35)

零到三阶非平衡态原点矩Hermite系数可写如下：

$\begin{split}&{\boldsymbol{a}}_1^{(0)} = \int {{f^{(1)}}{\rm{d}}{\boldsymbol{{{{\boldsymbol{\xi}} }} }}} = 0\\& {\boldsymbol{a}}_1^{(1)} = \int {{f^{(1)}}{\boldsymbol{{{{\boldsymbol{\xi}} }} }}{\rm{d}}{\boldsymbol{{{{\boldsymbol{\xi}} }} }}} = 0 \\& {\boldsymbol{a}}_1^{(2)} = \int {{f^{(1)}}({\boldsymbol{{{{\boldsymbol{\xi}} }} {{{\boldsymbol{\xi}} }} }} - {\boldsymbol{\delta }}){\rm{d}}{\boldsymbol{{{{\boldsymbol{\xi}} }} }}} \approx { \Pi ^{(1)} }\\& {\boldsymbol{a}}_1^{(3)} = \int {{f^{(1)}}({\boldsymbol{{{{\boldsymbol{\xi}} }} {{{\boldsymbol{\xi}} }} {{{\boldsymbol{\xi}} }} }} - 3{\boldsymbol{\delta {{{\boldsymbol{\xi}} }} }}){\rm{d}}{\boldsymbol{{{{\boldsymbol{\xi}} }} }}} \approx {Q^{(1)}}\\\end{split} $

(36)

其中二阶和三阶非平衡态系数可通过Chapman-Enskog展开^[73]得到：

$\begin{split}\Pi _{\alpha \beta }^{(1)} =& - \tau \rho \theta \Big({\partial _{\alpha} }{u_\beta } + {\partial _\beta }{u_{\alpha} } - \frac{2}{D}{\delta _{{\alpha} \beta }}{\partial _\gamma }{u_\gamma }\Big)\\ Q_{{\alpha} \beta \gamma }^{(1)} =& {u_{\alpha} }a_{1,\beta \gamma }^{(2)} + {u_\beta }a_{1,{\alpha} \gamma }^{(2)} + {u_\gamma }a_{1,{\alpha} \beta }^{(2)} -\\& \tau \rho \theta \Big({\delta _{{\alpha} \beta }}{\partial _\gamma }\theta + {\delta _{\beta \gamma }}{\partial _{\alpha} }\theta + {\delta _{{\alpha} \gamma }}{\partial _\beta }\theta \Big) \\\end{split}$

(37)

如果非平衡态分布函数用中心矩空间Hermite基进行展开^{[49, 55]}，则有：

$ {f^{({{1}})}}({\boldsymbol{x}},{\boldsymbol{{{{\boldsymbol{\xi}} }} }},t) \approx {f^{(1),N}}({\boldsymbol{x}},{\boldsymbol{{{{\boldsymbol{\xi}} }} }},t) = \omega ({\boldsymbol{c}})\sum\limits_{n = 0}^N {\frac{1}{{n!}}{\boldsymbol{b}}_1^{(n)}:{{\boldsymbol{H}}^{(n)}}({\boldsymbol{c}})} $

(38)

中心矩空间非平衡态分布函数的Hermite系数可计算如下：

$ {\boldsymbol{b}}_1^{(n)} = \int {{f^{(1)}}{{\boldsymbol{H}}^{(n)}}({\boldsymbol{c}}){\rm{d}}{\boldsymbol{c}}} $

(39)

利用中心矩Hermite基与原点矩Hermite基的数学关系，可将中心矩空间非平衡态的Hermite系数算出，下面直接给出零到四阶的显式表达式：

${\boldsymbol{b}}_1^{(0)} = 0 $

(40)

$ {\boldsymbol{b}}_1^{(1)} = 0 $

(41)

$ {\boldsymbol{b}}_1^{(2)} = {\boldsymbol{a}}_1^{(2)} $

(42)

$ {\boldsymbol{b}}_1^{(3)} = {\boldsymbol{a}}_1^{(3)} - 3{\boldsymbol{a}}_1^{(2)}{\boldsymbol{u}} $

(43)

$ {\boldsymbol{b}}_1^{(4)} = {\boldsymbol{a}}_1^{(4)} - 4{\boldsymbol{ua}}_1^{(3)} + 6{{\boldsymbol{u}}^2}{\boldsymbol{a}}_1^{(2)} $

(44)

亦可将原点矩非平衡Hermite系数用中心矩非平衡Hermite系数表示，显式写出零到四阶有：

$\begin{split}&{\boldsymbol{a}}_1^{(0)} = 0\\& {\boldsymbol{a}}_1^{(1)} = 0\\& {\boldsymbol{a}}_1^{(2)} = {\boldsymbol{b}}_1^{(2)}\\& {\boldsymbol{a}}_1^{(3)} = {\boldsymbol{b}}_1^{(3)} + 3{\boldsymbol{ub}}_1^{(2)}\\& {\boldsymbol{a}}_1^{(4)} = {\boldsymbol{b}}_1^{(4)} + 4{\boldsymbol{ub}}_1^{(3)} - 6{{\boldsymbol{u}}^2}{\boldsymbol{b}}_1^{(2)} \\\end{split}$

(45)

更进一步地，如果将非平衡态分布函数在温度标度的中心矩空间进行展开，则有：

$ {f^{({{1}})}}({\boldsymbol{x}},{\boldsymbol{{{{\boldsymbol{\xi}} }} }},t) \approx {f^{(1),N}}({\boldsymbol{x}},{\boldsymbol{{{{\boldsymbol{\xi}} }} }},t) = \omega ({\boldsymbol{\eta }})\sum\limits_{n = 0}^N {\frac{1}{{n!}}{\boldsymbol{d}}_1^{(n)}:{{\boldsymbol{H}}^{(n)}}({\boldsymbol{\eta }})} $

(46)

相应地，对应各阶Hermite系数可计算如下：

$ {\boldsymbol{d}}_1^{(n)} = \int {{f^{(1)}}{{\boldsymbol{H}}^{(n)}}({\boldsymbol{\eta }})\frac{1}{{{{\Big(\sqrt \theta \Big)}^D}}}{\rm{d}}{\boldsymbol{{{{\boldsymbol{\xi}} }} }}} $

(47)

利用温度标度中心矩Hermite基与原点矩Hermite基的数学关系，我们可以计算上式，显式地写出零到四阶：

$\begin{split}&{\boldsymbol{d}}_1^{(0)} = 0\\& {\boldsymbol{d}}_1^{(1)} = 0\\& {\boldsymbol{d}}_1^{(2)} = {\bigg(\frac{1}{{\sqrt \theta }}\bigg)^{D + 2}}{\boldsymbol{a}}_1^{(2)}\\& {\boldsymbol{d}}_1^{(3)} = {\bigg(\frac{1}{{\sqrt \theta }}\bigg)^{D + 3}}\Big({\boldsymbol{a}}_1^{(3)} - 3{\boldsymbol{ua}}_1^{(2)}\Big)\\& {\boldsymbol{d}}_1^{(4)} = {\bigg(\frac{1}{{\sqrt \theta }}\bigg)^{D + 4}}\Big[{\boldsymbol{a}}_1^{(4)} - 4{\boldsymbol{ua}}_1^{(3)} + 6({{\boldsymbol{u}}^2} - (\theta - 1){\boldsymbol{\delta }}){\boldsymbol{a}}_1^{(2)}\Big] \\\end{split}$

(48)

同时也可以用温度标度中心矩Hermite系数将原点矩Hermite系数表示为：

$\begin{split}&{\boldsymbol{a}}_1^{(0)} = 0 \\& {\boldsymbol{a}}_1^{(1)} = 0\\& {\boldsymbol{a}}_1^{(2)} = {\Big(\sqrt \theta \Big)^{D + 2}}{\boldsymbol{d}}_1^{(2)}\\& {\boldsymbol{a}}_1^{(3)} = {\Big(\sqrt \theta \Big)^{D + 3}}{\boldsymbol{d}}_1^{(3)} + 3{\boldsymbol{u}}{\Big(\sqrt \theta \Big)^{D + 2}}{\boldsymbol{d}}_1^{(2)}\\& {\boldsymbol{a}}_1^{(4)} = {\Big(\sqrt \theta \Big)^{D + 4}}{\boldsymbol{d}}_1^{(4)} + 4{\boldsymbol{u}}{\Big(\sqrt \theta \Big)^{D + 3}}{\boldsymbol{d}}_1^{(3)} -\\&\qquad\;\;\;\; 6\Big[{{\boldsymbol{u}}^2} + (\theta - 1){\boldsymbol{\delta }}\Big]{\Big(\sqrt \theta \Big)^{D + 2}}{\boldsymbol{d}}_1^{(2)}\\\end{split} $

(49)

1.4.3 力项的Hermite展开

式（1）中的力项为连续形式，下面对其在前N阶Hermite多项式构成Hilbert空间中进行讨论。力项的形式为 $ F({\boldsymbol{{{{\boldsymbol{\xi}} }} }}) = - {\boldsymbol{g}} \cdot {\nabla _{\boldsymbol{{{{\boldsymbol{\xi}} }} }}}{f_N} $ ，截断到N阶的粒子分布函数对粒子速度的导数可推导如下，其中利用了Hermite基的Rodrigues表达式（4）：

$ \begin{split}{\nabla _{\boldsymbol{{{{\boldsymbol{\xi}} }} }}}{f_N} \approx & \sum\limits_{n = 0}^N {\frac{1}{{n!}}{{\boldsymbol{a}}^{(n)}}:{\nabla _{\boldsymbol{{{{\boldsymbol{\xi}} }} }}}\Big[\omega ({\boldsymbol{{{{\boldsymbol{\xi}} }} }}){{\boldsymbol{H}}^{(n)}}({\boldsymbol{{{{\boldsymbol{\xi}} }} }})\Big]} \\=& - \omega ({\boldsymbol{{{{\boldsymbol{\xi}} }} }})\sum\limits_{n = 0}^N {\frac{1}{{n!}}{{\boldsymbol{a}}^{(n)}}:{{\boldsymbol{H}}^{(n + 1)}}({\boldsymbol{{{{\boldsymbol{\xi}} }} }})} \\ =& - \omega ({\boldsymbol{{{{\boldsymbol{\xi}} }} }})\sum\limits_{n = 1}^{N + 1} {\frac{1}{{n!}}\Big(n{{\boldsymbol{a}}^{(n - 1)}}\Big):{{\boldsymbol{H}}^{(n)}}({\boldsymbol{{{{\boldsymbol{\xi}} }} }})}\end{split} $

(50)

因此，力项的截断形式Hermite展开可写为：

$ F({\boldsymbol{{{{\boldsymbol{\xi}} }} }}) \approx \omega ({\boldsymbol{{{{\boldsymbol{\xi}} }} }})\sum\limits_{n = 1}^{N + 1} {\frac{1}{{n!}}\Big(n{\boldsymbol{g}}{{\boldsymbol{a}}^{(n - 1)}}\Big):{{\boldsymbol{H}}^{(n)}}({\boldsymbol{{{{\boldsymbol{\xi}} }} }})} $

(51)

其中 $ n{\boldsymbol{g}}{{\boldsymbol{a}}^{(n - 1)}} $ 为 $ {{\boldsymbol{a}}^{(n - 1)}} $ 和 $ {\boldsymbol{g}} $ 的并矢，亦为对称张量，如三阶形式为 $ 3{\boldsymbol{g}}{{\boldsymbol{a}}^{(2)}} = {g_{\alpha} }{a_{\beta \gamma }} + {g_\beta }{a_{\gamma {\alpha} }} + {g_\gamma }{a_{{\alpha} \beta }} $ 。力项在原点矩、中心矩和温度标度中心矩空间的展开可以分别写出：

$\begin{split}&F({\boldsymbol{{{{\boldsymbol{\xi}} }} }}) \approx \omega ({\boldsymbol{{{{\boldsymbol{\xi}} }} }})\sum\limits_{n = 1}^{N + 1} {\frac{1}{{n!}}n{\boldsymbol{g}}\Big({\boldsymbol{a}}_0^{(n - 1)} + {\boldsymbol{a}}_1^{(n - 1)}\Big):{{\boldsymbol{H}}^{(n)}}({\boldsymbol{{{{\boldsymbol{\xi}} }} }})}\\& F({\boldsymbol{{{{\boldsymbol{\xi}} }} }}) \approx \omega ({\boldsymbol{c}})\sum\limits_{n = 1}^{N + 1} {\frac{1}{{n!}}n{\boldsymbol{g}}\Big({\boldsymbol{b}}_0^{(n - 1)} + {\boldsymbol{b}}_1^{(n - 1)}\Big):{{\boldsymbol{H}}^{(n)}}({\boldsymbol{c}})}\\& F({\boldsymbol{{{{\boldsymbol{\xi}} }} }}) \approx \frac{{\omega ({\boldsymbol{\eta }})}}{{\sqrt \theta }}\sum\limits_{n = 1}^{N + 1} {\frac{1}{{n!}}n{\boldsymbol{g}}\Big({\boldsymbol{d}}_0^{(n - 1)} + {\boldsymbol{d}}_1^{(n - 1)}\Big):{{\boldsymbol{H}}^{(n)}}({\boldsymbol{\eta }})}\\\end{split} $

(52)

其中 $ {\nabla _{\boldsymbol{{{{\boldsymbol{\xi}} }} }}} = {\nabla _{\boldsymbol{c}}} $ 和 ${\nabla _{\boldsymbol{{{{\boldsymbol{\xi}} }} }}} = \Big({1 \mathord{\left/ {\vphantom {1 {\sqrt \theta }}} \right. } {\sqrt \theta }}\Big)^D{\nabla _{\boldsymbol{\eta }}}$ 关系式得到了应用。如果令N = 3（后文中会讨论非平衡态展开会比平衡态展开低一阶，且高阶项贡献很小），力项的截断显式展开表达式为：

$ \qquad\begin{split}F({\boldsymbol{{{{\boldsymbol{\xi}} }} }}) \approx& \omega ({\boldsymbol{{{{\boldsymbol{\xi}} }} }})\rho \bigg\{ \underbrace {{\boldsymbol{g}} \cdot {\boldsymbol{{{{\boldsymbol{\xi}} }} }}}_{1{\rm{st}}} + \underbrace {({\boldsymbol{g}} \cdot {\boldsymbol{{{{\boldsymbol{\xi}} }} }})({\boldsymbol{{{{\boldsymbol{\xi}} }} }} \cdot {\boldsymbol{u}}) - {\boldsymbol{g}} \cdot {\boldsymbol{u}}}_{2{\rm{nd}}} + \\&\underbrace {\frac{1}{{6\rho }}3{\boldsymbol{g}}\Big\{\rho \big[{{\boldsymbol{u}}^2} + (\theta - 1){\boldsymbol{\delta }}\big] + {\boldsymbol{a}}_1^{(2)}\Big\}:{{\boldsymbol{H}}^{(3)}}({\boldsymbol{{{{\boldsymbol{\xi}} }} }})}_{3{\rm{rd}}} +\\& \underbrace {\frac{1}{{24\rho }}4{\boldsymbol{g}}\Big({\boldsymbol{a}}_0^{(3)} + {\boldsymbol{a}}_1^{(3)}\Big):{{\boldsymbol{H}}^{(4)}}({\boldsymbol{{{{\boldsymbol{\xi}} }} }})}_{4{\rm{th}}}\bigg\} \end{split}$

(53)

$\qquad\begin{split} F({\boldsymbol{{{{\boldsymbol{\xi}} }} }}) \approx& \omega ({\boldsymbol{c}})\rho \bigg\{ \underbrace {{\boldsymbol{g}} \cdot {\boldsymbol{c}}}_{1{\rm{st}}} + \underbrace 0_{2{\rm{nd}}} + \\&\underbrace {\frac{1}{{6\rho }}3{\boldsymbol{g}}\Big[\rho (\theta - 1){\boldsymbol{\delta }} + {\boldsymbol{a}}_1^{(2)}\Big]:{{\boldsymbol{H}}^{(3)}}({\boldsymbol{c}})}_{3{\rm{rd}}} +\\& \underbrace {\frac{1}{{24\rho }}4{\boldsymbol{g}}\Big({\boldsymbol{a}}_1^{(3)} - 3{\boldsymbol{ua}}_1^{(2)}\Big):{{\boldsymbol{H}}^{(4)}}({\boldsymbol{c}})}_{4{\rm{th}}}\bigg\} \end{split}$

(54)

$\qquad \begin{split}F({\boldsymbol{{{{\boldsymbol{\xi}} }} }}) \approx& \frac{{\omega ({\boldsymbol{\eta }})}}{\Big({\sqrt \theta }\Big)^D}\rho \Bigg\{ \underbrace {{\boldsymbol{g}} \cdot {{\bigg(\frac{1}{{\sqrt \theta }}\bigg)}^D}{\boldsymbol{\eta }}}_{1{\rm{st}}} +\\& \underbrace 0_{2{\rm{nd}}} + \underbrace {\frac{1}{{6\rho }}3{\boldsymbol{g}}{{\bigg(\frac{1}{{\sqrt \theta }}\bigg)}^{D + 2}}{\boldsymbol{a}}_1^{(2)}:{{\boldsymbol{H}}^{(3)}}({\boldsymbol{\eta }})}_{3{\rm{rd}}} +\\& \underbrace {\frac{1}{{24\rho }}4{\boldsymbol{g}}{{\bigg(\frac{1}{{\sqrt \theta }}\bigg)}^{D + 3}}\Big({\boldsymbol{a}}_1^{(3)} - 3{\boldsymbol{ua}}_1^{(2)}\Big):{{\boldsymbol{H}}^{(4)}}({\boldsymbol{\eta }})}_{4{\rm{th}}}\Bigg\} \end{split}$

(55)

其中式（53）为Martys等推导的形式^[74]。分布函数的非平衡态Hermite矩在力项中的贡献十分微弱，因此在实际的连续流数值计算中常可忽略。He等^[75]的力项处理即为假设式（50）中分布函数由平衡态代替，得到了广泛应用。类似于分布函数，速度相空间力的离散形式为 $ {F_i} = {w_i}{{F({{\boldsymbol{{{{\boldsymbol{\xi}} }} }}_i})} \mathord{\left/ {\vphantom {{F({{\boldsymbol{{{{\boldsymbol{\xi}} }} }}_i})} {\omega ({{\boldsymbol{{{{\boldsymbol{\xi}} }} }}_i})}}} \right. } {\omega ({{\boldsymbol{{{{\boldsymbol{\xi}} }} }}_i})}} $ 。

令 $ {\boldsymbol{F}} = \rho {\boldsymbol{g}} $ ，力项展开的原点矩Hermite系数分别为：

$\begin{split}&{\boldsymbol{a}}_F^{(1)} = {\boldsymbol{F}}\\& {\boldsymbol{a}}_F^{(2)} = 2{\boldsymbol{Fu}}\\& {\boldsymbol{a}}_F^{(3)} = 3{\boldsymbol{F}}\big[{\boldsymbol{uu}} + (\theta - 1){\boldsymbol{\delta }}\big]\\& {\boldsymbol{a}}_F^{(4)} = 4{\boldsymbol{F}}\big[{\boldsymbol{uuu}} + 3(\theta - 1){\boldsymbol{u\delta }}\big]\\ \end{split}$

(56)

力项展开的中心矩Hermite系数为：

${\boldsymbol{b}}_F^{(1)} = {\boldsymbol{F}} $

(57)

${\boldsymbol{b}}_F^{(2)} = 0 $

(58)

${\boldsymbol{b}}_F^{(3)} = 3{\boldsymbol{F\delta }}(\theta - 1) $

(59)

$ {\boldsymbol{b}}_F^{(4)} = 0 $

(60)

力项展开的温度标度中心矩Hermite系数为：

$ \begin{split}&{\boldsymbol{d}}_F^{(1)} = \frac{{\boldsymbol{F}}}{{{{\Big(\sqrt \theta \Big)}^{2D}}}}\\& {\boldsymbol{d}}_F^{(2)} = 0 \\& {\boldsymbol{d}}_F^{(3)} = 0\\& {\boldsymbol{d}}_F^{(4)} = 0 \end{split}$

(61)

点评：由式（56）～式（61）可以看出，如果不考虑非平衡态粒子分布函数的影响，力项在各空间展开的Hermite系数与平衡态Hermite系数保持一致；尤其是温度标度中心矩空间中我们可以看到只有一阶的Hermite系数不等于零，力项只在动量方程层面上有影响。

如果将式（56）~式（61）只保留到二阶，那么我们可以得到如下的力项表达式：

$ F({\boldsymbol{{{{\boldsymbol{\xi}} }} }}) \approx \omega ({\boldsymbol{{{{\boldsymbol{\xi}} }} }})\rho \big[{\boldsymbol{g}} \cdot {\boldsymbol{{{{\boldsymbol{\xi}} }} }} + ({\boldsymbol{g}} \cdot {\boldsymbol{{{{\boldsymbol{\xi}} }} }})({\boldsymbol{{{{\boldsymbol{\xi}} }} }} \cdot {\boldsymbol{u}}) - {\boldsymbol{g}} \cdot {\boldsymbol{u}}\big] $

(62)

$F({\boldsymbol{{{{\boldsymbol{\xi}} }} }}) \approx \omega ({\boldsymbol{{{{\boldsymbol{\xi}} }} }} - {\boldsymbol{u}})\rho {\boldsymbol{g}} \cdot ({\boldsymbol{{{{\boldsymbol{\xi}} }} }} - {\boldsymbol{u}}) = {\boldsymbol{g}} \cdot ({\boldsymbol{{{{\boldsymbol{\xi}} }} }} - {\boldsymbol{u}}){f^{{\rm{eq}}}} = - {\boldsymbol{g}} \cdot {\nabla _{{\boldsymbol{{{{\boldsymbol{\xi}} }} }} - {\boldsymbol{u}}}}{f^{{\rm{eq}}}} $

(63)

$F({\boldsymbol{{{{\boldsymbol{\xi}} }} }}) \approx \omega ({\boldsymbol{\eta }}){\bigg(\frac{1}{{\sqrt \theta }}\bigg)^{2D}}\rho {\boldsymbol{g}} \cdot {\boldsymbol{\eta }} = {\boldsymbol{g}} \cdot{\bigg(\frac{1}{{\sqrt \theta }}\bigg)^D} {\boldsymbol{\eta }}{f^{{\rm{eq}}}} = - {\boldsymbol{g}} \cdot {\nabla _{\boldsymbol{\xi }}}{f^{{\rm{eq}}}} $

(64)

上面的推导中用到了Maxwell-Boltzmann平衡态分布函数，且推导式（63）中用到了等温假设。式（62）如果在速度空间进行离散且考虑离散格子效应，即为Guo等^[76]2002年推导的力项格式（具体离散在后文中阐述）；式（63）即为He等^[75]1998年设计的力项格式；式（64）与式（63）等价。

1.4.4 碰撞项的Hermite展开

关于碰撞项的展开，介绍两种投影方式。第一种为基于原点矩、中心矩和温度标度中心矩空间的Hermite基展开，即 ${{\boldsymbol{H}}^{(n)}}({\boldsymbol{{{{\boldsymbol{\xi}} }} }})、{{\boldsymbol{H}}^{(n)}}({\boldsymbol{c}})、{{\boldsymbol{H}}^{(n)}}({\boldsymbol{\eta }})$ ；第二种为基于原点矩、中心矩和温度标度中心矩空间的单项式幂次基展开，即 ${{\boldsymbol{{{{\boldsymbol{\xi}} }} }}^n}、{{\boldsymbol{c}}^n}、{{\boldsymbol{\eta }}^n}$ 。第二种在Chen 等^[46]的工作和Shan^[55]的工作中被采纳。为了方便比较和厘清各种模型之间的关系，我们对两种方式的展开进行论述。

基于Hermite基展开碰撞项，可把Boltzmann-BGK方程的右端项展开如下：

$ \varOmega = \omega ({\boldsymbol{\varsigma }})\sum\limits_{n = {{2}}}^N {\frac{1}{{n!}}{\boldsymbol{p}}_\varOmega ^{(n)}:{{\boldsymbol{H}}^{(n)}}({\boldsymbol{\varsigma }})} $

(65)

其中 ${\boldsymbol{p}}_\varOmega ^{(n)} = {\boldsymbol{a}}_\varOmega ^{(n)}、{\boldsymbol{b}}_\varOmega ^{(n)}、{\boldsymbol{d}}_\varOmega ^{(n)}$ ，分别对应三个不同空间中的碰撞项Hermite系数。为方便分析，非平衡态在三种不同空间中展开，我们再次写出：

$ {f^{(1)}} = \omega ({\boldsymbol{\varsigma }})\sum\limits_{n = {{2}}}^N {\frac{1}{{n!}}{\boldsymbol{q}}_1^{(n)}:{{\boldsymbol{H}}^{(n)}}({\boldsymbol{\varsigma }})} $

(66)

其中， ${\boldsymbol{q}}_1^{(n)} = {\boldsymbol{a}}_1^{(n)}、{\boldsymbol{b}}_1^{(n)}、{\boldsymbol{d}}_1^{(n)}$ ，而自变量可选为 ${\boldsymbol{\varsigma }} = {\boldsymbol{{{{\boldsymbol{\xi}} }} }}、{\boldsymbol{c}}、{\boldsymbol{\eta }}$ 任意一个。不同空间中，碰撞项的Hermite系数与非平衡态的Hermite系数之间的转换关系相同，这里直接沿用。正如Chen等^[46]指出，多松弛碰撞项必须满足某些约束条件才能保证输运系数的伽利略不变性，非平衡分布函数的矩松弛必须与原始的碰撞项的松弛在讨论的阶数上分别等价，即：

$ \int \varOmega {{\boldsymbol{\varsigma }}^n}{\rm{d}}{\boldsymbol{\varsigma }} = - \frac{1}{{{\tau _n}}}\int {{f^{(1)}}{{\boldsymbol{\varsigma }}^n}{\rm{d}}{\boldsymbol{\varsigma }}} $

(67)

其中多弛豫时间被应用，不是原始的Boltzmann-BGK方程中的单松弛因子，而是在每阶矩应用了不同的松弛因子，可视为其拓展模型。在我们的工作中，上式的约束条件被拓展到非平衡态分布函数的Hermite矩松弛，这必须与原始的碰撞项的松弛在讨论的阶数上分别等价，即：

$ \int \varOmega {{\boldsymbol{H}}^{(n)}}({\boldsymbol{\varsigma }}){\rm{d}}{\boldsymbol{\varsigma }} = - \frac{1}{{{\tau _n}}}\int {{f^{(1)}}{{\boldsymbol{H}}^{(n)}}({\boldsymbol{\varsigma }}){\rm{d}}{\boldsymbol{\varsigma }}} $

(68)

首先我们讨论约束条件式（67）。 $ {{\boldsymbol{\varsigma }}^n} $ 与 $ {{\boldsymbol{H}}^{(n)}}({\boldsymbol{\varsigma }}) $ 前几阶的关系可由式（21）得到，任意阶数学关系可参考文献[55]，写出零到四阶如下：

$\begin{split}&1 = {{\boldsymbol{H}}^{(0)}}({\boldsymbol{\varsigma }})\\& {\boldsymbol{\varsigma }} = {{\boldsymbol{H}}^{(1)}}({\boldsymbol{\varsigma }})\\& {\boldsymbol{\varsigma \varsigma }} = {{\boldsymbol{H}}^{(2)}}({\boldsymbol{\varsigma }}) + {\boldsymbol{\delta }}{{\boldsymbol{H}}^{(0)}}({\boldsymbol{\varsigma }})\\& {\boldsymbol{\varsigma \varsigma \varsigma }} = {{\boldsymbol{H}}^{(3)}}({\boldsymbol{\varsigma }}) + 3{\boldsymbol{\delta }}{{\boldsymbol{H}}^{(1)}}({\boldsymbol{\varsigma }})\\& {\boldsymbol{\varsigma \varsigma \varsigma \varsigma }} = {{\boldsymbol{H}}^{(4)}}({\boldsymbol{\varsigma }}) + 6{{\boldsymbol{H}}^{(2)}}({\boldsymbol{\varsigma }}){\boldsymbol{\delta }} + 3{\boldsymbol{\delta \delta }}{{\boldsymbol{H}}^{(0)}}({\boldsymbol{\varsigma }}) \\\end{split}$

(69)

将上式代入式（67），并将碰撞项和非平衡态分布函数的Hermite展开式也代入，并利用Hermite多项式的正交特性，可得到下面关系式：

$ \begin{split}&{\boldsymbol{p}}_\varOmega ^{(2)} = - {\omega _2}{\boldsymbol{q}}_1^{(2)}\\& {\boldsymbol{p}}_\varOmega ^{(3)} = - {\omega _3}{\boldsymbol{q}}_1^{(3)}\\& {\boldsymbol{p}}_\varOmega ^{(4)} + 6{\boldsymbol{\delta p}}_\varOmega ^{(2)} = - {\omega _4}\Big({\boldsymbol{q}}_1^{(4)} + 6{\boldsymbol{\delta q}}_1^{(2)}\Big)\\\end{split} $

(70)

如果令 $ {\boldsymbol{p}}_\varOmega ^{(n)} = {\boldsymbol{b}}_\varOmega ^{(n)} $ 和 $ {\boldsymbol{q}}_1^{(n)} = {\boldsymbol{b}}_1^{(n)} $ ，即在中心矩空间进行展开，那么利用非平衡态的中心矩Hermite系数与原点矩Hermite系数之间的关系式（40）～式（44）、式（45），代入上式，则有：

$\begin{split}& {\boldsymbol{a}}_\varOmega ^{(2)} = - {\omega _2}{\boldsymbol{b}}_1^{(2)} = - {\omega _2}{\boldsymbol{a}}_1^{(2)}\\& {\boldsymbol{a}}_\varOmega ^{(3)} - 3{\boldsymbol{ua}}_\omega ^{(2)} = - {\omega _3}{\boldsymbol{b}}_1^{(3)} = - {\omega _3}\Big({\boldsymbol{a}}_1^{(3)} - 3{\boldsymbol{ua}}_1^{(2)}\Big)\\& {\boldsymbol{a}}_\varOmega ^{(4)} - 4{\boldsymbol{ua}}_\varOmega ^{(3)} + 6({\boldsymbol{uu}} + {\boldsymbol{\delta }}){\boldsymbol{a}}_\varOmega ^{(2)} = - {\omega _4}\Big({\boldsymbol{b}}_1^{(4)} + 6{\boldsymbol{\delta b}}_1^{(2)}\Big) =\\& \qquad- {\omega _4}\Big[{\boldsymbol{a}}_1^{(4)} - 4{\boldsymbol{ua}}_1^{(3)} + 6({\boldsymbol{uu}} + {\boldsymbol{\delta }}){\boldsymbol{a}}_1^{(2)}\Big] \\\end{split}$

(71)

整理式（71）后的碰撞项为：

$ {\boldsymbol{a}}_\varOmega ^{(2)} = - {\omega _2}{\boldsymbol{a}}_1^{(2)} $

(72)

$ {\boldsymbol{a}}_\varOmega ^{(3)} = - {\omega _3}{\boldsymbol{a}}_1^{(3)} + 3({\omega _3} - {\omega _2}){\boldsymbol{ua}}_1^{(2)} $

(73)

$\begin{split} {\boldsymbol{a}}_\varOmega ^{(4)} =& - {\omega _4}{\boldsymbol{a}}_1^{(4)} + 4({\omega _4} - {\omega _3}){\boldsymbol{ua}}_1^{(3)} -\\& 6({\omega _4} + {\omega _2} - 2{\omega _3}){\boldsymbol{uua}}_1^{(2)} - 6({\omega _4} - {\omega _2}){\boldsymbol{\delta a}}_1^{(2)} \end{split}$

(74)

如果令 $ {\boldsymbol{p}}_\varOmega ^{(n)} = {\boldsymbol{d}}_\varOmega ^{(n)} $ 和 $ {\boldsymbol{q}}_1^{(n)} = {\boldsymbol{d}}_1^{(n)} $ ，进行上述类似的代数计算，发现最后得到的碰撞项Hermite系数与上述方程一模一样。这说明在该约束条件下，用不同的单项式幂次基对推导的碰撞项Hermite系数并无影响。碰撞项式（72）～式（74）的二阶和三阶与Chen等^[46]推导的碰撞形式等价，在Shan等^[47]于2007年提出的多松弛模型的基础上做了三阶项松弛的修正，即式（73）右端的第二项。Shan^[55]在2019年将Chen等的工作拓展到任意阶，并截断到四阶，即式（74）。

下面我们讨论利用约束条件式（68）得到碰撞项Hermite系数形式。将碰撞项和非平衡态分布函数的Hermite展开式代入式（68），并利用Hermite多项式的正交特性可以得到：

$ \begin{split}&{\boldsymbol{p}}_\varOmega ^{(2)} = - {\omega _2}{\boldsymbol{q}}_1^{(2)}\\& {\boldsymbol{p}}_\varOmega ^{(3)} = - {\omega _3}{\boldsymbol{q}}_1^{(3)}\\& {\boldsymbol{p}}_\varOmega ^{(4)} = - {\omega _4}{\boldsymbol{q}}_1^{(4)}\\\end{split} $

(75)

如果令 $ {\boldsymbol{p}}_\varOmega ^{(n)} = {\boldsymbol{a}}_\varOmega ^{(n)} $ 和 $ {\boldsymbol{q}}_1^{(n)} = {\boldsymbol{a}}_1^{(n)} $ ，则可得到原点矩空间的碰撞形式如下：

$\begin{split}&{\boldsymbol{a}}_\varOmega ^{(2)} = - {\omega _2}{\boldsymbol{a}}_1^{(2)}\\& {\boldsymbol{a}}_\varOmega ^{(3)} = - {\omega _3}{\boldsymbol{a}}_1^{(3)} \\& {\boldsymbol{a}}_\varOmega ^{(4)} = - {\omega _4}{\boldsymbol{a}}_1^{(4)}\\ \end{split}$

(76)

如果令 $ {\boldsymbol{p}}_\varOmega ^{(n)} = {\boldsymbol{b}}_\varOmega ^{(n)} $ 和 $ {\boldsymbol{q}}_1^{(n)} = {\boldsymbol{b}}_1^{(n)} $ ，并利用中心矩和原点矩Hermite系数之间的转换关系，可得到中心矩空间的碰撞形式：

$\begin{split}&{\boldsymbol{a}}_\varOmega ^{(2)} = {\boldsymbol{b}}_\varOmega ^{(2)} = - {\omega _2}{\boldsymbol{b}}_1^{(2)}\\& {\boldsymbol{a}}_\varOmega ^{(3)} - 3{\boldsymbol{ua}}_\varOmega ^{(2)} = {\boldsymbol{b}}_\varOmega ^{(3)} = - {\omega _3}{\boldsymbol{b}}_1^{(3)} \\& {\boldsymbol{a}}_\varOmega ^{(4)} - 4{\boldsymbol{ua}}_\varOmega ^{(3)} + 6{{\boldsymbol{u}}^2}{\boldsymbol{a}}_\varOmega ^{(2)} = {\boldsymbol{b}}_\varOmega ^{(4)} = - {\omega _4}{\boldsymbol{b}}_1^{(4)}\\ \end{split}$

(77)

或者进一步整理成下面形式：

${\boldsymbol{a}}_\varOmega ^{(2)} = - {\omega _2}{\boldsymbol{a}}_1^{(2)} $

(78)

${\boldsymbol{a}}_\varOmega ^{(3)} = - {\omega _3}{\boldsymbol{a}}_1^{(3)} + ({\omega _3} - {\omega _2})3{\boldsymbol{ua}}_1^{(2)} $

(79)

$\begin{split} {\boldsymbol{a}}_\varOmega ^{(4)} =& - {\omega _4}{\boldsymbol{a}}_1^{(4)} + ({\omega _4} - {\omega _3})4{\boldsymbol{ua}}_1^{(3)} - \\&({\omega _4} + {\omega _2} - 2{\omega _3})6{{\boldsymbol{u}}^2}{\boldsymbol{a}}_1^{(2)} \end{split}$

(80)

如果令 $ {\boldsymbol{p}}_\varOmega ^{(n)} = {\boldsymbol{d}}_\varOmega ^{(n)} $ 和 $ {\boldsymbol{q}}_1^{(n)} = {\boldsymbol{d}}_1^{(n)} $ ，并利用温度标度中心矩和原点矩Hermite系数之间的转换关系，可得到温度标度中心矩空间的碰撞形式：

$\begin{split}& {\boldsymbol{a}}_\varOmega ^{(2)} = {\Big(\sqrt \theta \Big)^{D + 2}}{\boldsymbol{d}}_\varOmega ^{(2)} = - {\omega _2}{\Big(\sqrt \theta \Big)^{D + 2}}{\boldsymbol{d}}_1^{(2)}\\& {\boldsymbol{a}}_\varOmega ^{(3)} - 3{\boldsymbol{ua}}_\varOmega ^{(2)} = {\Big(\sqrt \theta \Big)^{D + 3}}{\boldsymbol{d}}_\varOmega ^{(3)} = - {\omega _3}{\Big(\sqrt \theta \Big)^{D + 3}}{\boldsymbol{d}}_1^{(3)}\\& {\boldsymbol{a}}_\varOmega ^{(4)} - 4{\boldsymbol{ua}}_\varOmega ^{(3)} + 6{{\boldsymbol{u}}^2}{\boldsymbol{a}}_\varOmega ^{(2)} - 6(\theta - 1){\boldsymbol{\delta a}}_\varOmega ^{(2)} =\\&\qquad {\Big(\sqrt \theta \Big)^{D + 4}}{\boldsymbol{d}}_\varOmega ^{(4)} = - {\omega _4}{\Big(\sqrt \theta \Big)^{D + 4}}{\boldsymbol{d}}_1^{(4)} \\\end{split}$

(81)

或者整理成下面形式：

${\boldsymbol{a}}_\varOmega ^{(2)} = - {\omega _2}{\boldsymbol{a}}_1^{(2)} $

(82)

${\boldsymbol{a}}_\varOmega ^{(3)} = - {\omega _3}{\boldsymbol{a}}_1^{(3)} + ({\omega _3} - {\omega _2})3{\boldsymbol{ua}}_1^{(2)} $

(83)

$\begin{split} {\boldsymbol{a}}_\varOmega ^{(4)} =& - {\omega _4}{\boldsymbol{a}}_1^{(4)} + 4({\omega _4} - {\omega _3}){\boldsymbol{ua}}_1^{(3)} - 6({\omega _4} + {\omega _2} - \\&2{\omega _3}){{\boldsymbol{u}}^2}{\boldsymbol{a}}_1^{(2)} + 6({\omega _4} - {\omega _2})(\theta - 1){\boldsymbol{\delta a}}_1^{(2)}\end{split} $

(84)

式（82）～式（84）即为作者等在2019年提出的碰撞形式^[56]。比较碰撞形式式（78）～式（80）与式（82）～式（84），二阶三阶项无差别，四阶项中式（84）多了一项与温度偏差 $ (\theta - 1) $ 有关的修正项。

点评：从式（75）可以看出，不管是Hermite原点矩、中心矩还是温度标度中心矩，每阶都可以独立松弛。然而从本文1.4.2节对非平衡态Hermite系数在不同空间之间的数学关系的分析，我们可以发现，二阶矩在不同空间中是等价的，只要离散速度足够多，能保证二阶矩的积分精度，分子速度空间离散形式的非平衡态分布函数的二阶矩就会等价于连续形式非平衡态分布函数的二阶矩。Shan等^{[45, 71]}已证明对于二维至少需要17个离散速度，三维至少需要39个离散速度。只要离散速度充分，同时平衡态展开到三阶及以上，不论在哪个空间进行松弛，黏性输运系数的伽利略不变性都能得到严格满足^[55]。使用D2Q9或者D3Q27离散格子级联碰撞模型^[17]或者累积量碰撞模型^[27]都不能保证严格满足黏性输运系数的伽利略不变性。在不同的平移坐标系中，宏观速度是不一样的，且是有量纲的，在温度标度中心矩空间进行碰撞，各阶矩是与宏观速度无关且被局部声速无量纲化的量。从式（82）～式（84）即可看出，在温度标度中心矩空间进行松弛，松弛的量是真正的非平衡态量，在原点矩松弛对象中将守恒量（动量、内能）剔除，从而保证了每阶矩松弛的真正独立。结合式（76）、式（82）～式（84）我们可以发现，原点矩空间对每阶非平衡态Hermite矩进行松弛，实质上在三阶及以上矩松弛中与低阶矩发生了相互干扰，即所谓的串扰效应。这也说明原点矩空间中的各阶矩松弛并不真正的独立。因此，在碰撞算子中，松弛对象必须满足坐标的平移不变性。Li等^[77]最近的工作讨论了松弛对象必须满足坐标的旋转不变性。本综述中暂不作展开。

2 碰撞算子的正则化 2.1 格子玻尔兹曼方程的时空离散

包含力项的速度相空间离散形式的格子玻尔兹曼方程为：

$ \frac{{\partial {f_i}}}{{\partial t}} + {{{{\boldsymbol{\xi}} }} _i} \cdot \nabla {f_i} = {\varOmega _i} + {F_i} $

(85)

将上式沿特征线进行积分，右端积分采用梯形法计算，则有

$ \begin{split}&{f_i}({\boldsymbol{x}} + {{\boldsymbol{{{{\boldsymbol{\xi}} }} }}_i},t + 1) - {f_i}({\boldsymbol{x}},t) = \frac{1}{2}\big[{\varOmega _i}({\boldsymbol{x}} + {{\boldsymbol{{{{\boldsymbol{\xi}} }} }}_i},t + 1) + {\varOmega _i}({\boldsymbol{x}},t)\big] +\\&\qquad \frac{1}{2}\big[{F_i}({\boldsymbol{x}} + {{\boldsymbol{{{{\boldsymbol{\xi}} }} }}_i},t + 1) + {F_i}({\boldsymbol{x}},t)\big]\end{split} $

(86)

为简便起见，将时间步简化为 $ \Delta t = 1 $ 代入。上式为隐式方程，右端项包含下一时刻项，为消除上述方程的隐式特性，引入一个新定义变量：

$ \bar {{f_i}} = {f_i} - \frac{1}{2}{\varOmega _i} - \frac{1}{2}{F_i} $

(87)

将新定义变量代入式（86），可得到如下显式方程，右端项均为当前时间步的量：

$ \bar {{f_i}} ({\boldsymbol{x}} + {{\boldsymbol{{{{\boldsymbol{\xi}} }} }}_i},t + 1) = \bar {{f_p}} ({\boldsymbol{x}},t) = \bar {{f_i}} ({\boldsymbol{x}},t) + {\varOmega _i}({\boldsymbol{x}},t) + {F_i}({\boldsymbol{x}},t) $

(88)

其中 $\bar {{f_p}} ({\boldsymbol{x}},t)$ 代表碰撞后迁移前的粒子分布函数。将上式右端用Hermite展开，并乘以n阶Hermite基，即可得到右端项的Hermite系数，类似于式（6）操作。该操作本质上可视为粒子分布函数空间朝相应Hermite矩空间的一个转换，我们写出通用形式为：

$ \bar {\boldsymbol{q}} _p^{(n)} = {\bar {\boldsymbol{q}} ^{(n)}} - \frac{1}{{{\tau _n}}}\Big({{\boldsymbol{q}}^{(n)}} - {\boldsymbol{q}}_0^{(n)}\Big) + {\boldsymbol{q}}_F^{(n)} $

(89)

上式可为原点矩、中心矩或温度标度中心矩。

类似地，新定义变量转换到Hermite矩空间为：

$ {\bar {\boldsymbol{q}} ^{(n)}} = {{\boldsymbol{q}}^{(n)}} + \frac{1}{{2{\tau _n}}}\Big({{\boldsymbol{q}}^{(n)}} - {\boldsymbol{q}}_0^{(n)}\Big) - \frac{1}{2}{\boldsymbol{q}}_F^{(n)} $

(90)

因此，原始变量可用新定义变量表示为下式：

$ {{\boldsymbol{q}}^{(n)}} = \frac{{2{\tau _n}}}{{2{\tau _n} + 1}}{\bar {\boldsymbol{q}} ^{(n)}} + \frac{1}{{2{\tau _n} + 1}}{\boldsymbol{q}}_0^{(n)} + \frac{{{\tau _n}}}{{2{\tau _n} + 1}}{\boldsymbol{q}}_F^{(n)} $

(91)

将上式代入式（89），即为有关于新定义分布函数的矩空间松弛方程：

$ \bar {\boldsymbol{q}} _p^{(n)} = {\bar {\boldsymbol{q}} ^{(n)}} - \frac{1}{{{\lambda _n}}}\Big({\bar {\boldsymbol{q}} ^{(n)}} - {\boldsymbol{q}}_0^{(n)}\Big) + \Bigg(1 - \frac{1}{{2{\lambda _n}}}\Bigg){\boldsymbol{q}}_F^{(n)} $

(92)

其中上式中 $ {\lambda _n} = {\tau _n} + {1 \mathord{\left/ {\vphantom {1 2}} \right. } 2} $ 。新定义分布函数的非平衡态矩为：

$ \bar {\boldsymbol{q}} _1^{(n)} = {\bar {\boldsymbol{q}} ^{(n)}} - {\boldsymbol{q}}_0^{(n)} $

(93)

因此矩空间松弛方程亦可写为：

$ \bar {\boldsymbol{q}} _p^{(n)} = {\boldsymbol{q}}_0^{(n)} + \Bigg(1 - \frac{1}{{{\lambda _n}}}\Bigg)\bar {\boldsymbol{q}} _1^{(n)} + \Bigg(1 - \frac{1}{{2{\lambda _n}}}\Bigg){\boldsymbol{q}}_F^{(n)},\;\;\;\;n \geqslant 1 $

(94)

对应到原点矩、中心矩和温度标度中心矩则为：

$\begin{split}&\bar {\boldsymbol{a}} _p^{(n)} = {\boldsymbol{a}}_0^{(n)} + \Bigg(1 - \frac{1}{{{\lambda _n}}}\Bigg)\bar {\boldsymbol{a}} _1^{(n)} + \Bigg(1 - \frac{1}{{2{\lambda _n}}}\Bigg){\boldsymbol{a}}_F^{(n)},\;\;\;n \geqslant 1\\& \bar {\boldsymbol{b}} _p^{(n)} = {\boldsymbol{b}}_0^{(n)} + \Bigg(1 - \frac{1}{{{\lambda _n}}}\Bigg)\bar {\boldsymbol{b}} _1^{(n)} + \Bigg(1 - \frac{1}{{2{\lambda _n}}}\Bigg){\boldsymbol{b}}_F^{(n)},\;\;\;n \geqslant 1 \\& \bar {\boldsymbol{d}} _p^{(n)} = {\boldsymbol{d}}_0^{(n)} + \Bigg(1 - \frac{1}{{{\lambda _n}}}\Bigg)\bar {\boldsymbol{d}} _1^{(n)} + \Bigg(1 - \frac{1}{{2{\lambda _n}}}\Bigg){\boldsymbol{d}}_F^{(n)},\;\;\;n \geqslant 1 \\\end{split}$

(95)

特别注意地，当不包含力项时，碰撞方程从二阶开始，因为质量和动量守恒。但是对于新定义的分布函数，如果包含力项，则一阶非平衡态矩不为零，这在多相流^[62-63]和浸没边界法^[11]包含力项时需要注意。先对新定义的粒子分布函数所求得的宏观量与实际物理量之间的关系做一个讨论。

零阶矩为：

$ \sum {\bar {{f_i}} } = \sum {{f_i} - \frac{1}{2}{\varOmega _i} - \frac{1}{2}{F_i}} = \sum f $

(96)

可以得到实际密度为：

$ \rho = \bar \rho $

(97)

一阶矩为：

$ \sum {\bar {{f_i}} {{\boldsymbol{{{{\boldsymbol{\xi}} }} }}_i}} = \sum {\Bigg({f_i} - \frac{1}{2}{\varOmega _i} - \frac{1}{2}{F_i}\Bigg){{\boldsymbol{{{{\boldsymbol{\xi}} }} }}_i}} = \sum {\Bigg({f_i} - \frac{1}{2}{F_i}\Bigg){{\boldsymbol{{{{\boldsymbol{\xi}} }} }}_i}} $

(98)

可以得到实际速度为：

$ {\boldsymbol{u}} = \bar {\boldsymbol{u}} + \frac{{\boldsymbol{F}}}{{2\rho }} $

(99)

二阶矩的迹为：

$ \sum { {{{\bar f}_i}} {c_{i,{\alpha} }}{c_{i,{\alpha} }}} = \sum {\Bigg({f_i} - \frac{1}{2}{\varOmega _i} - \frac{1}{2}{F_i}\Bigg){c_{i,{\alpha} }}{c_{i,{\alpha} }}} = \sum {{f_i}{c_{i,{\alpha} }}{c_{i,{\alpha} }}} $

(100)

可以得到实际温度为：

$ \theta = \bar \theta $

(101)

在平衡态和力项中采用的均为实际物理量，即不带一拔标记的密度、速度和温度。而新定义的分布函数的一阶矩并不为零：

$ \begin{split}\bar {\boldsymbol{a}} _1^{(1)} = &\sum { {{\bar f}_i^{(1)}} {{\boldsymbol{{{{\boldsymbol{\xi}} }} }}_i}} = \sum {\Bigg(f_i^{(1)} - \frac{1}{2}{\varOmega _i} - \frac{1}{2}{F_i}\Bigg){{\boldsymbol{{{{\boldsymbol{\xi}} }} }}_i}} \\=& \sum { - \frac{1}{2}{F_i}{{\boldsymbol{{{{\boldsymbol{\xi}} }} }}_i}} = - \frac{1}{2}{\boldsymbol{F}} \end{split}$

(102)

类似地，对中心矩和温度标度中心矩一阶非平衡态可评估为：

$ \begin{split}\bar {\boldsymbol{b}} _1^{(1)} =& \sum { { {\bar f}_i^{(1)}} {{\boldsymbol{c}}_i}} = \sum {\Bigg(f_i^{(1)} - \frac{1}{2}{\varOmega _i} - \frac{1}{2}{F_i}\Bigg){{\boldsymbol{c}}_i}} \\ =&\sum { - \frac{1}{2}{F_i}{{\boldsymbol{c}}_i}} = - \frac{1}{2}{\boldsymbol{F}}\end{split} $

(103)

$\begin{split}\bar {\boldsymbol{d}} _1^{(1)} =& \frac{1}{{{{\Big(\sqrt \theta \Big)}^D}}}\sum { {{\bar f}_i^{(1)}} {{\boldsymbol{\eta }}_i}} = \frac{1}{{{{\Big(\sqrt \theta \Big)}^D}}}\sum {\Bigg(f_i^{(1)} - \frac{1}{2}{\varOmega _i} - \frac{1}{2}{F_i}\Bigg){{\boldsymbol{\eta }}_i}} \\ =&- \frac{1}{2}\frac{1}{{{{\Big(\sqrt \theta \Big)}^{D + 1}}}}{\boldsymbol{F}}\\[-15pt]\end{split} $

(104)

对于原点矩空间碰撞，平衡态部分、非平衡态部分以及力项均可以直接利用Hermite展开的截断阶形式重构粒子分布函数。而对于中心矩和温度标度中心矩空间的碰撞形式，由于离散速度依赖于当地速度和温度，我们选择将其转换到原点矩空间，再进行粒子分布函数重构。根据矩空间碰撞执行的空间不同以及截断阶的阶数，我们可以对已存在的正则化模型进行系统地分析和讨论。为方便起见，后文中非平衡态矩上面的一拔标记全部省略。

2.2 正则化碰撞模型分析与讨论 2.2.1 Ladd 正则化模型和Latt正则化模型

Ladd等在文献[39]（1994年）和文献[41]（2001年）中提出了正则化碰撞模型的最早形式。在他们的碰撞模型中，平衡态分布函数展开到二阶，非平衡态分布函数也仅保留到二阶。碰撞后的分布函数写为下面形式：

$ f_{p,i}^{{\rm{reg}}} = {\omega _i}\rho \Bigg[1 + {\boldsymbol{u}} \cdot {{\boldsymbol{{{{\boldsymbol{\xi}} }} }}_i} + \frac{{\Big({\boldsymbol{uu}} + {{{\boldsymbol{\Pi }}_p^{{\rm{neq}}}} \mathord{\left/ {\vphantom {{{\boldsymbol{\Pi }}_p^{{\rm{neq}}}} \rho }} \right. } \rho }\Big):\Big({{\boldsymbol{{{{\boldsymbol{\xi}} }} }}_i}{{\boldsymbol{{{{\boldsymbol{\xi}} }} }}_i} - c_s^2{\boldsymbol{I}}\Big)}}{2}\Bigg] $

(105)

注意：本文中速度量纲的量均被声速标度过了，因此声速不再出现在表达式中。上式的二阶非平衡态松弛为：

$ {\boldsymbol{\Pi }}_p^{{\rm{neq}}} = \Bigg(1 - \frac{1}{{{\lambda _{2,s}}}}\Bigg){\bar {\boldsymbol{\Pi }} ^{{\rm{neq}}}} + \Bigg(1 - \frac{1}{{{\lambda _{2,b}}}}\Bigg)\frac{1}{D}Tr({{\boldsymbol{\Pi }}^{{\rm{neq}}}}){\boldsymbol{I}} $

(106)

${\bar {\boldsymbol{\Pi }} ^{{\rm{neq}}}}$ 为二阶应力张量的无迹部分， $Tr({{\boldsymbol{\Pi }}^{{\rm{neq}}}})$ 为二阶应力张量的迹， $ {\lambda _{2,s}} $ 为剪切粘性弛豫时间， $ {\lambda _{2,b}} $ 为体积粘性弛豫时间。二阶应力张量由非平衡态粒子分布函数的二阶矩求得，如下式：

$ \Pi _{{\alpha} \beta }^{{\rm{neq}}} = \sum\limits_i {f_i^{(1)}{{{{\boldsymbol{\xi}} }} _{i,{\alpha} }}{{{{\boldsymbol{\xi}} }} _{i,\beta }}} $

(107)

如果上式中的离散速度为D2Q9、D3Q15、D3Q19、D3Q27等标准格子，则上式求得的二阶应力张量是有迹的，而非平衡态粒子分布函数的二阶积分是无迹的。即，如果精确计算，实际是没有体积粘性项的，之所以二阶应力张量中出现有迹部分，是因为离散速度不够多。体积粘性通常由多原子气体的内部自由度激发而产生，或者考虑稠密流体效应。Ladd的模型为等温模型，考虑到能量守恒，体积粘性模态的本征值（ $ {1 \mathord{\left/ {\vphantom {1 {{\lambda _{2,b}}}}} \right. } {{\lambda _{2,b}}}} $ ）应为零；但实际数值计算中不考虑能量守恒，该松弛因子自由可调。

Latt等^[42]在2005年也提出了正则化模型，他们的模型并未分离二阶应力张量的无迹部分和有迹部分的独立松弛，仅为

$ {\boldsymbol{\Pi }}_p^{{\rm{neq}}} = \Bigg(1 - \frac{1}{{{\lambda _2}}}\Bigg){{\boldsymbol{\Pi }}^{{\rm{neq}}}} $

(108)

Latt的模型在Ladd模型之后很久才提出，而且可视为Ladd模型的一个特例，即剪切和体积粘性松弛因子相等，本质上不能算一个新的模型。Ladd模型为原点矩空间的正则化碰撞模型，高于三阶的超出等温Navier-Stokes方程的非平衡态矩可视为松弛因子都设置为1，对碰撞后的粒子分布函数不产生贡献。

2.2.2 Zhang-Shan-Chen正则化模型

2006年，Zhang、Shan和Chen^[43]基于Hermite展开提出了高阶正则化模型，并用来模拟高Kn数非连续流动。他们认为平衡态分布函数投影到有限阶Hermite基构成的Hilbert空间，非平衡态部分也应该投影到该Hilbert空间上，并引入力项。根据本文2.1节的讨论，在考虑外力项后，非平衡态Hermite系数的一阶项并不为零，因此正则化的碰撞后粒子分布函数应该写为：

$ \begin{split}f_{p,i}^{{\rm{reg}}} = f_i^{(0)} + {\omega _i}\sum\limits_{i = 1}^3 {\frac{1}{{n!}}\Big({\boldsymbol{a}}_1^{(n)} + {\boldsymbol{a}}_\varOmega ^{(n)}\Big)} :{{\boldsymbol{H}}^{(n)}}({{\boldsymbol{{{{\boldsymbol{\xi}} }} }}_i}) + {F_i}\end{split} $

(109)

$\begin{split}f_{p,i}^{{\rm{reg}}} =& {\omega _i}\sum\limits_{n = 0}^3 {\frac{1}{{n!}}{\boldsymbol{a}}_0^{(n)}:{\boldsymbol{H}}({{\boldsymbol{{{{\boldsymbol{\xi}} }} }}_i})} + {\omega _i}\sum\limits_{n = 1}^3 {\Bigg(1 - \frac{1}{{{\lambda _n}}}\Bigg)\frac{1}{{n!}}{\boldsymbol{a}}_1^{(n)}:{\boldsymbol{H}}({{\boldsymbol{{{{\boldsymbol{\xi}} }} }}_i})} + \\& {\omega _i}\sum\limits_{n = 1}^3 {\Bigg(1 - \frac{1}{{2{\lambda _n}}}\Bigg)\frac{1}{{n!}}{\boldsymbol{a}}_F^{(n)}:{\boldsymbol{H}}({{\boldsymbol{{{{\boldsymbol{\xi}} }} }}_i})} \end{split}$

(110)

对于大于等于四阶的非平衡态矩（即 $ n \geqslant 4 $ ）全部截断，我们也可视为高阶非平衡态矩松弛因子全部设为 $ {\lambda _n} = 1 $ 。式（110）中一阶非平衡态松弛和力项的Hermite系数为：

$ \Bigg(1 - \frac{1}{{{\lambda _1}}}\Bigg)\Bigg( - \frac{{\boldsymbol{F}}}{2}\Bigg) + \Bigg(1 - \frac{1}{{2{\lambda _1}}}\Bigg){\boldsymbol{F}} = \frac{{\boldsymbol{F}}}{2} $

(111)

由上式可见，一阶松弛因子并无实际意义。

然而，在原始的Zhang-Shan-Chen模型中，式（109）中的重构是从二阶开始，且力项前面缺少系数 $ 1 - {1 \mathord{\left/ {\vphantom {1 {(2{\lambda _n})}}} \right. } {(2{\lambda _n})}} $ ，因此该模型理论上并不严谨。

2.2.3 Mattila-Philippi-Hegele正则化模型

Mattila等^[49]提出了一个基于中心矩空间展开的高阶格子正则化模型来处理可压缩传热流动，我们利用本文1.4.4节的有关碰撞项和1.4.2节的非平衡态系数在中心矩空间和原点矩空间分别展开的数学关系，分析该正则化模型。该模型中不包含力项，因此我们也省略力项。为分析方便，我们把非平衡态Hermite系数和碰撞项合并到一起，如式（95）所示，在式（77）中引入新的碰撞形式，包含非平衡的中心矩Hermite系数，具体如下：

$ {\boldsymbol{a}}_{\varOmega ,s}^{(2)} = {\boldsymbol{b}}_{\varOmega ,s}^{(2)} = (1 - {\omega _2}){\boldsymbol{b}}_1^{(2)} $

(112)

$ {\boldsymbol{a}}_{\varOmega ,s}^{(3)} - 3{\boldsymbol{ua}}_{\varOmega ,s}^{(2)} = {\boldsymbol{b}}_{\varOmega ,s}^{(3)} = (1 - {\omega _3}){\boldsymbol{b}}_1^{(3)} $

(113)

$ {\boldsymbol{a}}_{\varOmega ,s}^{(4)} - 4{\boldsymbol{ua}}_{\varOmega ,s}^{(3)} + 6{{\boldsymbol{u}}^2}{\boldsymbol{a}}_{\varOmega ,s}^{(2)} = {\boldsymbol{b}}_{\varOmega ,s}^{(4)} = (1 - {\omega _4}){\boldsymbol{b}}_1^{(4)} = 0 $

(114)

其中四阶及以上非平衡态松弛因子设为1，处理可压缩传热流动。即认为超出二阶应力张量和三阶热流矢量的高阶动理学中心矩快速松弛，碰撞后的粒子分布函数直接为平衡态。利用Hermite中心矩和原点矩之间的关系可整理上式为：

$ {\boldsymbol{a}}_{\varOmega ,s}^{(2)} = {\boldsymbol{a}}_1^{(2)} - {\omega _2}{\boldsymbol{a}}_1^{(2)} = {\boldsymbol{a}}_1^{(2)} + {\boldsymbol{a}}_\varOmega ^{(2)} $

(115)

${\boldsymbol{a}}_{\varOmega ,s}^{(3)} = {\boldsymbol{a}}_1^{(3)} + \Big[ - {\omega _3}{\boldsymbol{a}}_1^{(3)} + ({\omega _3} - {\omega _2})3{\boldsymbol{ua}}_1^{(2)}\Big] = {\boldsymbol{a}}_1^{(3)} + {\boldsymbol{a}}_\varOmega ^{(3)} $

(116)

$ \begin{split}{\boldsymbol{a}}_{\varOmega ,s}^{(4)} =& 4{\boldsymbol{ua}}_{\varOmega ,s}^{(3)} - 6{{\boldsymbol{u}}^2}{\boldsymbol{a}}_{\varOmega ,s}^{(2)} = \Big(4{\boldsymbol{ua}}_1^{(3)} - 6{{\boldsymbol{u}}^2}{\boldsymbol{a}}_1^{(2)}\Big) +\\& 4{\boldsymbol{u}}\Big[ - {\omega _3}{\boldsymbol{a}}_1^{(3)} + ({\omega _3} - {\omega _2})3{\boldsymbol{ua}}_1^{(2)}\Big] - 6{{\boldsymbol{u}}^2}\Big( - {\omega _2}{\boldsymbol{a}}_1^{(2)}\Big) \\=& \Big[4{\boldsymbol{ua}}_1^{(3)} - 6{{\boldsymbol{u}}^2}{\boldsymbol{a}}_1^{(2)}\Big] + 4{\boldsymbol{ua}}_\varOmega ^{(3)} - 6{{\boldsymbol{u}}^2}{\boldsymbol{a}}_\varOmega ^{(2)}\\[-12pt]\end{split} $

(117)

由上式（115）、式（116）可见，新的碰撞形式的二阶和三阶即为原碰撞项与非平衡态项之和，而由式（44）可知，去除四阶非平衡态中心矩 $ {\boldsymbol{b}}_1^{(4)} $ 后的四阶非平衡态原点矩即为 $ {\boldsymbol{a}}_1^{(4)} = 4{\boldsymbol{ua}}_1^{(3)} - 6{{\boldsymbol{u}}^2}{\boldsymbol{a}}_1^{(2)} $ ，恰为式（117）第一项。而将式（80）中的 $ {\boldsymbol{a}}_1^{(4)} $ 由去除中心矩后的形式 $ {\boldsymbol{a}}_1^{(4)} = 4{\boldsymbol{ua}}_1^{(3)} - 6{{\boldsymbol{u}}^2}{\boldsymbol{a}}_1^{(2)} $ 代入后，即为式（117）形式。因此，在本文建立的理论框架下，所谓的正则化过程即为在碰撞项式（114）中令四阶及以上的中心矩松弛因子为1，令高阶动理学矩快速到达平衡态。

如果是处理等温流动，可将三阶及以上的中心矩松弛因子都置为1，即有：

$\begin{split}&{\boldsymbol{a}}_{\varOmega ,s}^{(2)} = {\boldsymbol{b}}_{\varOmega ,s}^{(2)} = (1 - {\omega _2}){\boldsymbol{b}}_1^{(2)}\\& {\boldsymbol{a}}_{\varOmega ,s}^{(3)} - 3{\boldsymbol{ua}}_{\varOmega ,s}^{(2)} = {\boldsymbol{b}}_{\varOmega ,s}^{(3)} = (1 - {\omega _3}){\boldsymbol{b}}_1^{(3)} = 0\\& {\boldsymbol{a}}_{\varOmega ,s}^{(4)} - 4{\boldsymbol{ua}}_{\varOmega ,s}^{(3)} + 6{{\boldsymbol{u}}^2}{\boldsymbol{a}}_{\varOmega ,s}^{(2)} = {\boldsymbol{b}}_{\varOmega ,s}^{(4)} = (1 - {\omega _4}){\boldsymbol{b}}_1^{(4)} = 0\\\end{split} $

(118)

经过代数计算，整理上式为：

${\boldsymbol{a}}_{\varOmega ,s}^{(2)} = {\boldsymbol{a}}_1^{(2)} - {\omega _2}{\boldsymbol{a}}_1^{(2)} = {\boldsymbol{a}}_1^{(2)} + {\boldsymbol{a}}_\varOmega ^{(2)} $

(119)

$ {\boldsymbol{a}}_{\varOmega ,s}^{(3)} = 3{\boldsymbol{ua}}_{\varOmega ,s}^{(2)} = 3{\boldsymbol{ua}}_1^{(2)} - {\omega _2}(3{\boldsymbol{ua}}_1^{(2)}) $

(120)

${\boldsymbol{a}}_{\varOmega ,s}^{(4)} = 4{\boldsymbol{ua}}_{\varOmega ,s}^{(3)} - 6{{\boldsymbol{u}}^2}{\boldsymbol{a}}_{\varOmega ,s}^{(2)} = 6{{\boldsymbol{u}}^2}{\boldsymbol{a}}_{\varOmega ,s}^{(2)} = 6{{\boldsymbol{u}}^2}{\boldsymbol{a}}_1^{(2)} + 6{{\boldsymbol{u}}^2}\Big( - {\omega _2}{\boldsymbol{a}}_1^{(2)}\Big) $

(121)

K. Mattila等^[49]认为，在模拟可压缩传热流动时，平衡态保留到四阶矩，非平衡态也保留到四阶矩；但非平衡态部分 $ {\boldsymbol{a}}_1^{(4)} = {\boldsymbol{b}}_1^{(4)} + 4{\boldsymbol{ua}}_1^{(3)} - 6{{\boldsymbol{u}}^2}{\boldsymbol{a}}_1^{(2)} $ 中的中心矩是所谓的耗散部分，不能依赖于五阶平衡态矩，需要被滤除，即四阶非平衡态矩只保留所谓的对流部分： $ {\boldsymbol{a}}_1^{(4)} = 4{\boldsymbol{ua}}_1^{(3)} - 6{{\boldsymbol{u}}^2}{\boldsymbol{a}}_1^{(2)} $ 。写成下标形式为：

$\begin{split} {\boldsymbol{b}}_{1,a}^{(4)} =& 4{\boldsymbol{ua}}_1^{(3)} - 6{{\boldsymbol{u}}^2}{\boldsymbol{a}}_1^{(2)} =\Big({u_{\alpha} }a_{1,\beta \gamma \delta }^{(3)} + {u_\beta }a_{1,{\alpha} \gamma \delta }^{(3)} +\\& {u_\gamma }a_{1,{\alpha} \beta \delta }^{(3)} + {u_\delta }a_{1,{\alpha} \beta \gamma }^{(3)}\Big) - \Big({u_{\alpha} }{u_\beta }a_{1,\gamma \delta }^{(2)} + {u_{\alpha} }{u_\gamma }a_{1,\beta \delta }^{(2)} + \\&{u_{\alpha} }{u_\delta }a_{1,\beta \gamma }^{(2)} + {u_\beta }{u_\gamma }a_{1,{\alpha} \delta }^{(2)} + {u_\beta }{u_\delta }a_{1,{\alpha} \gamma }^{(2)} + {u_\gamma }{u_\delta }a_{1,{\alpha} \beta }^{(2)}\Big) \end{split}$

(122)

Mattila等的处理本质上与式（114）等价，在中心矩空间进行正则化，将四阶及以上的松弛因子置为1，也可视为滤除四阶及以上非平衡态中心矩。Zhang-Shan-Chen正则化模型则是在原点矩空间将四阶及以上松弛因子置为1， $ {\boldsymbol{a}}_1^{(4)} $ 的贡献全部滤掉；而在中心矩空间进行正则化， $ {\boldsymbol{a}}_1^{(4)} = 4{\boldsymbol{ua}}_1^{(3)} - 6{{\boldsymbol{u}}^2}{\boldsymbol{a}}_1^{(2)} $ ，四阶非平衡态原点矩的一部分被保留，该部分由低阶非平衡态矩和速度可递归算出。

在等温层级上，三阶及以上的非平衡态中心矩松弛因子均置为1，即均被滤除：

$ {\boldsymbol{b}}_{1,a}^{(3)} = {\boldsymbol{a}}_1^{(3)} - {\boldsymbol{b}}_1^{(3)} = 3{\boldsymbol{ua}}_1^{(2)} = {u_{\alpha} }a_{1,\beta \gamma }^{(2)} + {u_\beta }a_{1,{\alpha} \gamma }^{(2)} + {u_\gamma }a_{1,{\alpha} \beta }^{(2)} $

(123)

$\begin{split} {\boldsymbol{b}}_{1,a}^{(4)} =& 6{{\boldsymbol{u}}^2}{\boldsymbol{a}}_1^{(2)} = {u_{\alpha} }{u_\beta }a_{1,\gamma \delta }^{(2)} + {u_{\alpha} }{u_\gamma }a_{1,\beta \delta }^{(2)} + {u_{\alpha} }{u_\delta }a_{1,\beta \gamma }^{(2)} +\\& {u_\beta }{u_\gamma }a_{1,{\alpha} \delta }^{(2)} + {u_\beta }{u_\delta }a_{1,{\alpha} \gamma }^{(2)} + {u_\gamma }{u_\delta }a_{1,{\alpha} \beta }^{(2)} \end{split}$

(124)

式（120）、式（121）两式的右端第一项与上两式等价，而右端第二项为对上两式的松弛。另外一个值得注意的是，在Mattila的正则化模型中，单松弛碰撞形式 $ {\boldsymbol{a}}_\varOmega ^{(n)} = - {{{\boldsymbol{a}}_1^{(n)}} \mathord{\left/ {\vphantom {{{\boldsymbol{a}}_1^{(n)}} \lambda }} \right. } \lambda } $ 被采纳，只能处理普朗特数为1的流动，即使采用多松弛形式 $ {\boldsymbol{a}}_\varOmega ^{(n)} = - {{{\boldsymbol{a}}_1^{(n)}} \mathord{\left/ {\vphantom {{{\boldsymbol{a}}_1^{(n)}} {{\lambda _n}}}} \right. } {{\lambda _n}}} $ ，也不能实现热传导系数的伽利略不变性。这在Chen等的工作^[46]和Shan等的工作^[55]中已经通过数值测试证明过，从式（116）中就可以看出必须减除串扰项 $({\omega _3} - {\omega _2}) 3{\boldsymbol{ua}}_1^{(2)}$ 。

2.2.4 Malaspinas正则化模型

Malaspinas ^[48]在2015年提出了一个等温层级的递归正则化模型。递归包括平衡态Hermite系数的递归以及非平衡态Hermite系数的递归。平衡态Hermite系数的递归关系式在只考虑等温Navier-Stokes方程时很简单，可由数学归纳法得到：

$ {\boldsymbol{a}}_0^{(n)} = {\boldsymbol{a}}_0^{(n - 1)}{\boldsymbol{u}} {，} {\boldsymbol{a}}_0^{(0)} = \rho $

(125)

非平衡态Hermite系数的递归关系根据平衡态Hermite系数的递归关系、零阶CE展开得到的欧拉方程和一阶CE展开得到的非平衡态Hermite系数与平衡态Hermite系数时空导数之间的关系，推导得到：

$ a_{1,{{\alpha} _1} \ldots {{\alpha} _n}}^{(n)} = a_{1,{{\alpha} _1} \ldots {{\alpha} _{n - 1}}}^{(n - 1)}{u_{{{\alpha} _n}}} + \Big[{u_{{{\alpha} _1}}} \ldots {u_{{{\alpha} _{n - 2}}}}a_{1,{{\alpha} _{n - 1}}{{\alpha} _n}}^{(2)} + perm({{\alpha} _n})\Big],n \geqslant 3 $

(126)

其中， $ perm({{\alpha} _n}) $ 代表除 $ {{\alpha} _n} $ 外的索引的全排列，上式代表了n阶非平衡态Hermite系数可由n-1阶、二阶非平衡态Hermite系数和局部速度递归得到。

对于二维问题，他在Maxwell-Boltzmann平衡态的二阶泰勒展开的基础上，给出了部分三阶和四阶项，系数由递归关系给出。对于D2Q9格子，他延续了D'Humières^[15-16]和Dellar^[78]等的思想，认为有多少个离散速度就应该有多少个基，下面展开式除了前面零到二阶的六个基。另外还包含三阶和四阶的三个基：

$ \begin{split} f_i^{(0)} =& {\omega _i}\rho \Big[1 + {{\boldsymbol{{{{\boldsymbol{\xi}} }} }}_i} \cdot {\boldsymbol{u}} + \frac{1}{2}H_i^{(2)}:{\boldsymbol{uu}} + \\& \frac{1}{2}\Big(H_{i,xxy}^{(3)}u_x^2{u_y} + H_{i,xyy}^{(3)}u_y^2{u_x}\Big) + \frac{1}{4}H_{i,xxyy}^{(4)}u_x^2u_y^2\Big] \end{split} $

(127)

相应地，非平衡态展开为：

$ \begin{split}f_i^{(1)} =& {\omega _i}\bigg[\frac{1}{2}H_i^{(2)}:{\boldsymbol{a}}_1^{(2)} + \frac{1}{2}\Big(H_{i,xxy}^{(3)}a_{1,xxy}^{(3)} + H_{i,xyy}^{(3)}a_{1,xyy}^{(3)}\Big) +\\& \frac{1}{4}H_{i,xxyy}^{(4)}a_{1,xxyy}^{(4)}\bigg]\\[-12pt] \end{split}$

(128)

展开基跟平衡态展开基对应，二阶非平衡态Hermite系数通过非平衡态粒子分布函数的二阶矩求得，三阶和四阶非平衡态Hermite系数则递归求得。

相应地，三维（D3Q27）的平衡态展开零到二阶跟二维一样，但高阶一直展开到六阶，一共27个Hermite基，见式（129）：

$ \begin{split}f_i^{(0)} =& {\omega _i}\rho \bigg\{ 1 + {{\boldsymbol{{{{\boldsymbol{\xi}} }} }}_i} \cdot {\boldsymbol{u}} + \frac{1}{2}H_i^{(2)}:{\boldsymbol{uu}} +\\& \frac{1}{2}\Big(H_{i,xxy}^{(3)}u_x^2{u_y} + H_{i,xxz}^{(3)}u_x^2{u_z} + H_{i,xyy}^{(3)}{u_x}u_y^2 + H_{i,xzz}^{(3)}{u_x}u_z^2 + H_{i,yzz}^{(3)}{u_y}u_z^2 + H_{i,yyz}^{(3)}u_y^2{u_z} + 2H_{i,xyz}^{(3)}{u_x}{u_y}{u_z}\Big) + \\& \frac{1}{4}\bigg[H_{i,xxyy}^{(4)}u_x^2u_y^2 + H_{i,xxzz}^{(4)}u_x^2u_z^2 + H_{i,yyzz}^{(4)}u_y^2u_z^2 + 2\Big(H_{i,xyzz}^{(4)}{u_x}{u_y}u_z^2 + H_{i,xyyz}^{(4)}{u_x}u_y^2{u_z} + H_{i,xxyz}^{(4)}u_x^2{u_y}{u_z}\Big)\bigg] +\\& \frac{1}{4}\Big(H_{i,xxyzz}^{(5)}u_x^2{u_y}u_z^2 + H_{i,xxyyz}^{(5)}u_x^2u_y^2{u_z} + H_{i,xyyzz}^{(5)}{u_x}u_y^2u_z^2\Big) + \frac{1}{8}H_{i,xxyyzz}^{(6)}u_x^2u_y^2u_z^2\bigg\}\end{split} $

(129)

非平衡态粒子分布函数采用一样的Hermite基展开，二阶以上非平衡态Hermite系数可递归得到，见式（130）：

$ \begin{split}f_i^{(1)} =& {\omega _i}\bigg\{ \frac{1}{2}H_i^{(2)}:{\boldsymbol{a}}_1^{(2)} +\\& \frac{1}{2}\Big(H_{i,xxy}^{(3)}a_{1,xxy}^{(3)} + H_{i,xxz}^{(3)}a_{1,xxz}^{(3)} + H_{i,xyy}^{(3)}a_{1,xyy}^{(3)} + H_{i,xzz}^{(3)}a_{1,xzz}^{(3)} + H_{i,yzz}^{(3)}a_{1,yzz}^{(3)} + H_{i,yyz}^{(3)}a_{1,yyz}^{(3)} + 2H_{i,xyz}^{(3)}a_{1,xyz}^{(3)}\Big) +\\& \frac{1}{4}\bigg[H_{i,xxyy}^{(4)}a_{1,xxyy}^{(4)} + H_{i,xxzz}^{(4)}a_{1,xxzz}^{(4)} + H_{i,yyzz}^{(4)}a_{1,yyzz}^{(4)} + 2\Big(H_{i,xyzz}^{(4)}a_{1,xyzz}^{(4)} + H_{i,xyyz}^{(4)}a_{1,xyyz}^{(4)} + H_{i,xxyz}^{(4)}a_{1,xxyz}^{(4)}\Big)\bigg] +\\& \frac{1}{4}\Big(H_{i,xxyzz}^{(5)}a_{1,xxyzz}^{(5)} + H_{i,xxyyz}^{(5)}a_{1,xxyyz}^{(5)} + H_{i,xyyzz}^{(5)}a_{1,xyyzz}^{(5)}\Big) + \frac{1}{8}H_{i,xxyyzz}^{(6)}a_{1,xxyyzz}^{(6)}\bigg\}\end{split} $

(130)

最后的正则化碰撞和迁移格式为：

$ {f_i}({\boldsymbol{x}} + {{\boldsymbol{{{{\boldsymbol{\xi}} }} }}_i},t + 1) = f_i^{(0)} + \Bigg(1 - \frac{1}{\lambda }\Bigg)f_i^{(1)} $

(131)

对于该模型，递归得到的三阶及以上非平衡态动理学矩的松弛与二阶应力张量矩的松弛因子保持一致，与Ladd^[39]和Latt^[42]的模型在平衡态构建和非平衡态重构、松弛上存在高阶差异。另外，重构平衡态和非平衡态粒子分布函数所用的基与级联碰撞模型^[19]中由一维张量积得到的基是一致的。等温层面上，非平衡态Hermite系数的递归关系得到的高阶非平衡态矩与K. Mattila等^[49]基于中心矩展开得到的正则化模型中的非平衡态Hermite系数是等价的，如三阶和四阶表达式（123）、式（124）与递归表达式（126）得到的系数是一致的；而更高阶项可将非平衡态Hermite原点矩和中心矩之间的关系表达式（45）展开到高阶，然后滤除三阶及以上的中心矩，即可得到四阶以上的原点矩非平衡态Hermite系数，与递归表达式（126）保持一致。可见，Malaspinas推导得到的非平衡态Hermite系数递归关系式也可以用中心矩展开和滤除二阶以上中心矩非平衡态Hermite系数得到。根据本文1.2节的CE展开，我们可以发现，二阶的非平衡Hermite系数的k阶量依赖于平衡态Hermite系数的k+2阶量，如 $ {\boldsymbol{a}}_1^{(2)} $ 将最高依赖于 $ {\boldsymbol{a}}_0^{(3)} $ ， $ {\boldsymbol{a}}_2^{(2)} $ 将最高依赖于 $ {\boldsymbol{a}}_0^{(4)} $ ，而 $ {\boldsymbol{a}}_3^{(2)} $ 则最高依赖于 $ {\boldsymbol{a}}_0^{(5)} $ ，即平衡态的Hermite展开需截断到五阶。因此，在等温N-S方程模拟中，通常平衡态Hermite展开只截断三阶，而不是无穷展开下去。这样利用非平衡态粒子分布函数的二阶矩求得的二阶非平衡态Hermite系数中，其实也隐藏着离散速度不够充分的幽灵模态。2018年，Jacob等^[51]又在Malasipinas递归正则化模型上做了一些改进，他们将二阶非平衡态系数 $ {\boldsymbol{a}}_1^{(2)} $ 用有限差分进行了计算，并与用非平衡态粒子分布函数的二阶矩求得的 $ {\boldsymbol{a}}_1^{(2)} $ 进行加权，在D3Q19格子上进行数值执行，得到所谓的混合递归正则化模型，并对其是否可控进行了探讨。该混合递归正则化模型属于引入了额外的超粘性，增强格式的数值稳定性。2019年，Chen ^[79]提出了另外一个正则化模型，认为他们的模型与Malaspinas的递归正则化模型是等价的，但目前数学等价关系尚不明朗。

2.2.5 Li-Shi-Shan正则化模型

Li等^[56]在2019年提出了一个基于温度标度中心矩空间展开的高阶格子正则化模型来处理可压缩传热流动。我们利用本文1.4.4节的有关碰撞项和1.4.2节的非平衡态系数在温度标度中心矩空间和原点矩空间分别展开的数学关系来分析该正则化模型。为了跟本文2.2.3节分析保持一致，我们把温度标度中心矩空间的非平衡态Hermite系数及其松弛进行合并，在式（81）的基础上定义新的碰撞形式：

$ {\boldsymbol{a}}_{\varOmega ,s}^{(2)} = {\Big(\sqrt \theta \Big)^{D + 2}}{\boldsymbol{d}}_{\varOmega ,s}^{(2)} = (1 - {\omega _2}){\Big(\sqrt \theta \Big)^{D + 2}}{\boldsymbol{d}}_1^{(2)} $

(132)

${\boldsymbol{a}}_{\varOmega ,s}^{(3)} - 3{\boldsymbol{ua}}_{\varOmega ,s}^{(2)} = {\Big(\sqrt \theta \Big)^{D + 3}}{\boldsymbol{d}}_{\varOmega ,s}^{(3)} = (1 - {\omega _3}){\Big(\sqrt \theta \Big)^{D + 3}}{\boldsymbol{d}}_1^{(3)} $

(133)

$\begin{split}&{\boldsymbol{a}}_{\varOmega ,s}^{(4)} - 4{\boldsymbol{ua}}_{\varOmega ,s}^{(3)} + 6{{\boldsymbol{u}}^2}{\boldsymbol{a}}_{\varOmega ,s}^{(2)} - 6(\theta - 1){\boldsymbol{\delta a}}_{\varOmega ,s}^{(2)} = \\&\qquad{\Big(\sqrt \theta \Big)^{D + 4}}{\boldsymbol{d}}_{\varOmega ,s}^{(4)} = (1 - {\omega _4}){\Big(\sqrt \theta \Big)^{D + 4}}{\boldsymbol{d}}_1^{(4)} = 0 \end{split}$

(134)

利用温度标度中心矩和原点矩空间非平衡态Hermite系数之间的关系式，可将上式整理为：

$\begin{split}{\boldsymbol{a}}_{\varOmega ,s}^{(2)} =& (1 - {\omega _2}){\boldsymbol{a}}_1^{(2)} \\ {\boldsymbol{a}}_{\varOmega ,s}^{(3)} = &{\boldsymbol{a}}_1^{(3)} + \Big[ - {\omega _3}{\boldsymbol{a}}_1^{(3)} + ({\omega _3} - {\omega _2})3{\boldsymbol{ua}}_1^{(2)}\Big] = {\boldsymbol{a}}_1^{(3)} + {\boldsymbol{a}}_\varOmega ^{(3)} \\ {\boldsymbol{a}}_{\varOmega ,s}^{(4)} =& 4{\boldsymbol{ua}}_{\varOmega ,s}^{(3)} - 6{{\boldsymbol{u}}^2}{\boldsymbol{a}}_{\varOmega ,s}^{(2)} + 6(\theta - 1){\boldsymbol{\delta a}}_{\varOmega ,s}^{(2)} = \Big(4{\boldsymbol{ua}}_1^{(3)} - 6{{\boldsymbol{u}}_1}{\boldsymbol{a}}_1^{(2)}\Big) +\\& 4{\boldsymbol{u}}\Big[ - {\omega _3}{\boldsymbol{a}}_1^{(3)} + ({\varOmega _3} - {\omega _2})3{\boldsymbol{ua}}_1^{(2)}\Big] - 6{{\boldsymbol{u}}_1}( - {\omega _2}{\boldsymbol{a}}_1^{(2)}) \\=& \bigg\{4{\boldsymbol{ua}}_1^{(3)} - 6\Big[{{\boldsymbol{u}}^2} + (\theta - 1){\boldsymbol{\delta }}\Big]{\boldsymbol{a}}_1^{(2)}\bigg\} + 4{\boldsymbol{ua}}_\varOmega ^{(3)} - 6{{\boldsymbol{u}}_1}{\boldsymbol{a}}_\varOmega ^{(2)} \end{split} $

(135)

其中 $ {{\boldsymbol{u}}_1} = {{\boldsymbol{u}}^2} + (\theta - 1){\boldsymbol{\delta }} $ 。从上式可以发现，二阶和三阶碰撞就是在碰撞形式（78）、式（79）上将原点矩Hermite非平衡态系数进行合并；而四阶矩碰撞中原点矩Hermite非平衡态系数则为滤除四阶温度标度中心矩Hermite系数 $ {\boldsymbol{d}}_1^{(4)} $ 之后的形式：

$ {\boldsymbol{a}}_1^{(4)} = 4{\boldsymbol{ua}}_1^{(3)} - 6\Big[{{\boldsymbol{u}}^2} + (\theta - 1){\boldsymbol{\delta }}\Big]{\boldsymbol{a}}_1^{(2)} $

(136)

如果将该形式代入式（80）中，进行代数整理，与式（134）一致。其最后的碰撞后粒子分布函数形式为：

$ f_i^p = f_i^{(0)} + {\omega _i}\sum\limits_{n = 2}^4 {\frac{1}{{n!}}{\boldsymbol{a}}_{\varOmega ,s}^{(n)}:{{\boldsymbol{H}}^{(n)}}({{\boldsymbol{{{{\boldsymbol{\xi}} }} }}_i})} $

(137)

因此，Li等提出的基于温度标度中心矩空间中Hermite展开的碰撞形式可视为将该空间中四阶及以上的非平衡态Hermite矩的松弛因子置为1，该空间中的高阶矩快速松弛到平衡态。比较Zhang-Shan-Chen在原点矩空间的正则化碰撞模型、Mattila-Philippi-Hegele在中心矩空间的正则化碰撞模型，以及作者等的Li-Shi-Shan在温度标度中心矩空间的正则化碰撞模型，三者的相同点是将四阶及以上的相应非平衡态高阶矩的松弛因子置为1，让这些高阶非平衡态动力学矩快速松弛到平衡态，而本质差异是进行碰撞的空间不同。这些模型的联系和差异，在本文的分析框架下变得十分明朗。另外，在温度标度中心矩空间的正则化模型，在包含力项后的形式，我们将在另外的论文中发表。

2.2.6 Coreixas等正则化模型

2017年，Coreixas等^[50]将Malaspinas提出的等温层级的递归正则化模型推广到了可压缩传热层级上，发展了一个高阶格子递归正则化模型。该正则化模型包含平衡态Hermite系数的递归关系式：

$ a_{0,{{\alpha} _1} \cdots {{\alpha} _n}}^{(n)} = {u_{{{\alpha} _n}}}a_{0,{{\alpha} _1} \cdots {{\alpha} _{n - 1}}}^{(n - 1)} + (\theta - 1)\sum\limits_{i = 1}^{n - 1} {{\delta _{{{\alpha} _i}{{\alpha} _n}}}a_{0,{\beta _i}}^{(n - 2)}} ,\;\;n \geqslant 2 $

(138)

其中， $ a_{0,{\beta _i}}^{(n)} = a_{0,{{\alpha} _1} \cdots {{\alpha} _{i - 1}}{{\alpha} _{i + 1}} \cdots {{\alpha} _n}}^{(n)} $ ， $ {\beta _i} $ 代表索引中没有下标 $ {{\alpha} _i} $ 。上式与Malaspinas提出的等温递归正则化模型相比，平衡态递归关系式多了与温度有关的第二项。非平衡态Hermite系数的递归关系式也是根据平衡态Hermite系数递归关系式、零阶和一阶CE展开推导得到，具体如下：

$\begin{split} a_{1,{{\alpha} _1} \cdots {{\alpha} _n}}^{(n)} =& {u_{{{\alpha} _n}}}a_{1,{{\alpha} _1} \cdots {{\alpha} _{n - 1}}}^{(n - 1)} + (\theta - 1)\sum\limits_{i = 1}^{n - 1} {{\delta _{{{\alpha} _i}{{\alpha} _n}}}a_{1,{\beta _i}}^{(n - 2)}} +\\& \frac{1}{\rho }\sum\limits_{i = 1}^{n - 1} {a_{0,{\beta _i}}^{(n - 2)}a_{1,{{\alpha} _i}{{\alpha} _n}}^{(2)}} + \frac{1}{\rho }\sum\limits_{i = 1}^{n - 1} {\sum\limits_{j > i}^{n - 1} {a_{0,{\beta _{ij}}}^{(n - 3)}} \cdot} \\&{\Big(a_{1,{{\alpha} _i}{{\alpha} _j}{{\alpha} _n}}^{(3)} - {u_{{{\alpha} _i}}}a_{1,{{\alpha} _j}{{\alpha} _n}}^{(2)} - {u_{{{\alpha} _j}}}a_{1,{{\alpha} _i}{{\alpha} _n}}^{(2)} - {u_{{{\alpha} _n}}}a_{1,{{\alpha} _i}{{\alpha} _j}}^{(2)}\Big)，} n \geqslant 4 \end{split}$

(139)

上述递归关系式十分复杂，其数学推导过程也极其复杂。我们对四阶非平衡态Hermite系数进行整理如下：

$\begin{split} a_{1,{\alpha} \beta \gamma \delta }^{(4)} =& {u_\delta }a_{1,{\alpha} \beta \gamma }^{(3)} + (\theta - 1)\sum\limits_{i = 1}^3 {{\delta _{{{\alpha} _i}\delta }}a_{1,{\beta _i}}^{(2)}} + \frac{1}{\rho }\sum\limits_{i = 1}^3 {a_{0,{\beta _i}}^{(2)}a_{1,{{\alpha} _i}\delta }^{(2)}} + \frac{1}{\rho }\Bigg[\sum\limits_{i = 1}^3 {\sum\limits_{j > i}^3 {a_{0,{\beta _{ij}}}^{(1)}\Big(a_{1,{{\alpha} _i}{{\alpha} _j}\delta }^{(3)} - {u_{{{\alpha} _i}}}a_{1,{{\alpha} _j}\delta }^{(2)} - {u_{{{\alpha} _j}}}a_{1,{{\alpha} _i}\delta }^{(2)} - {u_\delta }a_{1,{{\alpha} _i}{{\alpha} _j}}^{(2)}\Big)} } \Bigg] \\=& {u_\delta }a_{1,{\alpha} \beta \gamma }^{(3)} + (\theta - 1)\Big({\delta _{{\alpha} \delta }}a_{1,\beta \gamma }^{(2)} + {\delta _{\beta \delta }}a_{1,{\alpha} \gamma }^{(2)} + {\delta _{\gamma \delta }}a_{1,{\alpha} \beta }^{(2)}\Big) + \frac{1}{\rho }\Big(a_{0,\beta \gamma }^{(2)}a_{1,{\alpha} \delta }^{(2)} + a_{0,{\alpha} \gamma }^{(2)}a_{1,\beta \delta }^{(2)} + a_{0,{\alpha} \beta }^{(2)}a_{1,\gamma \delta }^{(2)}\Big) + \\& \frac{1}{\rho }\bigg[a_{0,\gamma }^{(1)}\Big(a_{1,{\alpha} \beta \delta }^{(3)} - {u_{\alpha} }a_{1,\beta \delta }^{(2)} - {u_\beta }a_{1,{\alpha} \delta }^{(2)} - {u_\delta }a_{1,{\alpha} \beta }^{(2)}\Big) + a_{0,\beta }^{(1)}\Big(a_{1,{\alpha} \gamma \delta }^{(3)} - {u_{\alpha} }a_{1,\gamma \delta }^{(2)} - {u_\gamma }a_{1,{\alpha} \delta }^{(2)} - {u_\delta }a_{1,{\alpha} \gamma }^{(2)}\Big) +\\& a_{0,{\alpha} }^{(1)}\Big(a_{1,\beta \gamma \delta }^{(3)} - {u_\beta }a_{1,\gamma \delta }^{(2)} - {u_\gamma }a_{1,\beta \delta }^{(2)} - {u_\delta }a_{1,\beta \gamma }^{(2)}\Big)\bigg] \\= &4{\boldsymbol{ua}}_1^{(3)} - 6\Big[{{\boldsymbol{u}}^2} - (\theta - 1){\boldsymbol{\delta }}\Big]{\boldsymbol{a}}_1^{(2)}\end{split} $

(140)

可以发现，上式与Li-Shi-Shan模型基于温度标度中心矩空间推导出来的四阶非平衡态矩式（135）是等价的。其实更高阶的非平衡态Hermite系数可通过类似于式（48）、式（49）的数学转换关系式将相应阶的 $ {\boldsymbol{d}}_1^{(n)} $ 滤除即可得到。从Mattila-Philippi-Hegele正则化碰撞模型与Malaspinas的递归正则化模型之间的联系，以及Li-Shi-Shan高阶正则化模型与Coreixas等的高阶递归正则化模型之间的联系，我们可以发现：等温层级的原点矩空间非平衡态Hermite系数之间的递归关系很容易由中心矩与原点矩空间Hermite系数之间数学转换关系得到；类似地，可压缩传热层级的原点矩空间非平衡态Hermite系数之间的递归关系很容易由温度标度中心矩与原点矩空间Hermite系数之间的数学转换关系得到。

在Coreixas等的正则化碰撞模型中，重构完非平衡态Hermite系数后，其碰撞项形式为：

$ {\boldsymbol{a}}_\varOmega ^{(n)} = - {{{\boldsymbol{a}}_1^{(n)}} \mathord{\left/ {\vphantom {{{\boldsymbol{a}}_1^{(n)}} {{\lambda _n}}}} \right. } {{\lambda _n}}} $

(141)

最后的碰撞后粒子分布函数重构表达式为：

$ f_i^p = f_i^{(0)} + {\omega _i}\sum\limits_{n = 2}^4 {\frac{1}{{n!}}\Bigg(1 - \frac{1}{{{\lambda _n}}}\Bigg){\boldsymbol{a}}_1^{(n)}:{{\boldsymbol{H}}^{(n)}}({{\boldsymbol{{{{\boldsymbol{\xi}} }} }}_i})} $

(142)

与Li-Shi-Shan正则化碰撞模型式（132）～式（134）相比，Coreixas等的高阶递归正则化模型即使使用了多松弛模型，也不能保证热传导系数的伽利略不变性；而Li-Shi-Shan正则化碰撞模型不仅可以得到非平衡态Hermite系数的递归形式，还能得到碰撞项的递归形式。这充分体现了对温度标度中心矩空间进行展开和正则化在数学和物理基本原理上的优势。

3 结　论

本文对格子玻尔兹曼领域的正则化碰撞模型进行了系统的严谨的理论回顾，通过采用不同自变量的Hermite基对平衡态、非平衡态、力项和碰撞项进行展开，将所有的理论模型在该理论体系下的定位以及它们之间的联系进行了分析。分析发现：

1） 每一种正则化碰撞模型都可以在Hermite展开的理论框架下得到。

2） Malaspinas 和Coreixas的递归正则化模型，非平衡态可以无限递归得到，对于超出所研究的动理学矩（二阶应力张量矩和三阶热流矢量矩），其弛豫时间不一定为1。前者由与剪切粘性相关的弛豫时间统一松弛，后者则对每一阶非平衡态原点矩进行独立松弛。高阶矩的松弛因子如何影响模型的数值稳定性和精度有待进一步探究。

3） 除2）中提到的递归正则化模型，其他的正则化模型等价于在原点矩空间、中心矩空间和温度标度中心矩空间，将高于所关注的非平衡态动理学矩的弛豫时间置为1。

4） 可模拟可压缩传热NSF方程的正则化模型中，只有Chen等^[46]的模型、Shan的模型^[55]、Li-Shi-Shan模型^[56]是在中心矩或温度标度中心矩空间进行松弛，满足热传导系数的伽利略不变性。

正则化多松弛碰撞模型的理论本质就是将离散速度表达不了的高阶动理学矩进行截断，或者等价于将高阶动理学矩的松弛因子置为1，以避免虚假模态对水动力学方程的负面影响；而高阶动理学矩的松弛因子不为1时，对数值稳定性的正面影响目前尚不明朗。展望未来，正则化模型相较于其他碰撞模型的优势还有待严格地理论证明，尤其是在湍流模拟、声学模拟、多相流以及稀薄气体流动等领域，还有很大的发展空间。