0 引　言

在过去的三十年间，LBM 已经被发展到诸多领域，除了作为可压、不可压Navier-Stokes （N-S）方程的替代工具，还可以用作其他偏微分方程的解法器。尤其是在不可压流动模拟方面，LBM 被广泛认为是一种简单高效的数值工具。标准LBM的碰撞迁移过程使得编程简单、易于并行、数值耗散低。然而，这样的求解方式也存在一些不足。例如，物理空间和离散速度空间的耦合导致标准LBM必须使用均匀笛卡尔网格，对复杂外形的几何适应性差，时间步长和空间步长必须一致，即Courant-Friedrichs-Lewy（CFL）数为1，标准的单松弛BGK（Bhatnagar-Gross-Krook）格式在高雷诺数情况下数值稳定性差。此外，尽管LBM的数值耗散低，然而它始终还是只有二阶空间精度。应用LBM进行湍流模拟往往需要大量的网格，计算效率低。提高标准LBM的精度是较为困难的，目前主要的策略是对速度空间高阶离散或者对平衡态函数高阶展开。这些不足给LBM应用到高速、复杂外形的工程问题造成了一定困难。

格子玻尔兹曼方程（LBE）可以看作是离散速度玻尔兹曼方程（DVBE）的一种差分格式。已证明，物理对称性和格子对称性可以分开处理^[1]，物理空间的离散可以独立于速度空间的离散，这样就可以将离散速度和计算网格进行解耦。在这一思想指导下建立起来的方法一般称为非标准格子玻尔兹曼方法（OLBM）。OLBM主要分为两大类，一类是在欧拉框架下采用传统的数值离散方法，比如利用有限差分^[2]、有限体积^[3]、有限元法^[4]直接求解DVBE（以下称为欧拉法）；另一类是在半拉格朗日（semi-Lagrangian）框架下借助插值方法实现粒子分布函数的演化（以下称为半拉格朗日法）^[5-7]。早期的OLBM大多只有二阶精度，并且相比标准LBM会引入更大的数值耗散。近十多年来，高精度格式开始引起学者们的广泛关注。这类格式拥有较小的数值耗散，在多尺度流动的模拟上具有明显优势，例如湍流、气动噪声。和低精度的格式相比，高精度格式在取得相同误差条件下需要的网格要少得多，因而提高了计算效率。利用高精度格式求解N-S方程数值模拟不可压流的算法一般较为复杂，主要因为N-S方程的对流项是非线性，黏性项包含二阶偏导数，以及压强和速度的耦合。最常用的投影法或者压强修正法都需要求解压强泊松方程，该方程不容易收敛。非线性的对流项对于通过插值基函数在网格单元内构造数值解高阶多项式的高精度方法来说，还可能会产生混淆误差。而DVBE可以看成一组线性偏微分方程，它的曲线坐标形式相比N-S方程也要简洁得多，所以用高精度格式求解DVBE要容易得多。事实上，相关工作早在21世纪初就有展开。Shi等^[8]采用间断伽辽金谱元法求解DVBE；Düster等^[9]也发展了一种间断伽辽金格子玻尔兹曼法。然而，由于这两种方法缺少网格收敛性研究，它们的收敛阶能否取得高阶并不清楚。后来，Hejranfar和Ezzatneshan^[10-11]提出了一种紧致有限差分格子玻尔兹曼法（CFDLBM）并做了收敛性分析，验证了该方法能够取得四阶精度。由于中心型的紧致有限差分格式没有数值耗散，为了确保计算的稳定，必须引入过滤器。他们将CFDLBM和标准LBM做了比较，当取得相似的误差时，CFDLBM的效率更高。为了避免过滤操作、简化计算过程，后续又有相关工作采用带有数值耗散的高阶差分格式。如利用五阶WENO格式^[12]以及五阶迎风型紧致有限差分格式^[13]求解DVBE。有限差分对于复杂外形的几何适应性较差，为了改进这一问题，Li^[14]将高精度谱差分方法（SD）和DVBE结合，提出了SDLBM。通过坐标变换形函数，SDLBM可以应用非结构网格。Li对SDLBM做了详细的收敛性分析，证明了在被压缩性误差影响之前，SDLBM可以获得高阶精度，并且算法简单、易于并行。

以上这些方法的时间离散都采用显式格式。当雷诺数较高时，DVBE碰撞项的刚性将严重限制时间步长，大大降低计算效率。目前应对这种问题的方法主要有两种，一种是将碰撞项隐式处理，直接求解DVBE（以下称为直接法），常用的方法有Guo和Zhao^[15]提出的基于半隐格式的显式方法和隐式-显式Runge-Kutta方法（IMEX）等^[16]；另一种是将计算分为两步（以下称为分步法），首先执行碰撞步，这和标准LBM一样，然后在迁移步采用高精度格式求解一个纯对流方程。通过这种方法学者们也提出了不同的高精度OLBM，如Min和Lee^[17]提出的谱元间断伽辽金格子玻尔兹曼法以及Zadehgol等^[18]提出的节点间断伽辽金格子玻尔兹曼法。然而，鲜有文献比较分步法和碰撞项隐式处理的直接法之间的数值表现。

另一种构造高精度OLBM的思路是在半拉格朗日框架下采用高阶插值^[19]。由于不需要求解对流方程，没有CFL条件的限制，这种方法能在大时间步长下保持稳定，并且每个时间步空间算子（插值）只演化一次，因此这种方法计算效率较高。然而事实上，半拉格朗日法误差和时间步长是成正比的。要想获得高精度，还是需要取较小的时间步长。此外，半拉格朗日法的离散误差相比分步法要更大^[20]。因此，这几种高精度OLBM的优劣依然值得进一步讨论。

本文将最近较为流行的通量重构（FR）格式引入OLBM，发展了一种高阶通量重构格子玻尔兹曼方法（FRLBM）。FR格式为多种单元内多自由度的高阶格式，如节点间断有限元（NDG）和SD，提供了一个统一的框架^[21]，且算法更加简单，计算效率更高^[22]。在FRLBM的框架下，本文比较了直接法和分步法的精度以及稳定性，并且对比了基于半隐格式的显式方法和一阶IMEX、二阶IMEX以及三阶IMEX对不同雷诺数下CFL数取值的影响。研究结果旨在进一步加深对高精度OLBM的认识，为发展更加高效高精度的OLBM提供参考。

1 FRLBM计算方法 1.1 控制方程

采用BGK近似的DVBE为：

$ \frac{{\partial {f_{\alpha} }}}{{\partial t}} + {{\boldsymbol{e}}_{\alpha} }\cdot\nabla {f_{\alpha} } = - \frac{{\left( {{f_{\alpha} } - f_{\alpha} ^{{\text{eq}}}} \right)}}{\tau } $

(1)

其中， $ {f_{\alpha} } $ 为粒子分布函数，下标 $ \alpha $ 表示离散的粒子速度方向； ${{\boldsymbol{e}}_{\alpha} }$ 为粒子速度矢量； $ \tau $ 为松弛时间；为了减少压缩性误差，采用He-Luo不可压模型^[23]，平衡态分布函数 $f_{\alpha} ^{{\rm{eq}}}$ 定义为：

$ f_{\alpha} ^{{\text{eq}}} = {\omega _{\alpha} }\left[ {p + {p_0}\left( {\frac{{{{\boldsymbol{e}}_{\alpha} }\cdot{\boldsymbol{u}}}}{{c_s^2}} + \frac{{{{\left( {{{\boldsymbol{e}}_{\alpha} }\cdot{\boldsymbol{u}}} \right)}^2}}}{{2c_s^4}} - \frac{{{\boldsymbol{u}}\cdot{\boldsymbol{u}}}}{{2c_s^2}}} \right)} \right] $

(2)

其中， ${\omega _{\alpha} }$ 为权系数； $ {p_0} $ 为参考压强，直接设为1；流体宏观压强 $ p $ 和宏观速度 ${\boldsymbol{u}}$ 可以由下式求得：

$ p = \sum\limits {{f_{\alpha} }} ,\quad {p_0}{\boldsymbol{u}} = \sum\limits {{{\boldsymbol{e}}_{\alpha} }} {f_{\alpha} } $

(3)

本文采用D2Q9离散速度模型，相应的离散速度和权系数为：

$ {{\boldsymbol{e}}_{\alpha} } = \left\{ {\begin{array}{*{20}{l}} {(0,0),}&{\alpha = 0} \\ {( \pm 1,0),(0, \pm 1),}&{\alpha = 1 ～ 4} \\ {( \pm 1, \pm 1),}&{\alpha = 5 ～ 8} \end{array}} \right. $

(4)

$ {\omega _{\alpha} } = \left\{ {\begin{array}{*{20}{l}} {4/9,}&{\alpha = 0} \\ {1/9,}&{\alpha = 1 ～ 4} \\ {1/36,}&{\alpha = 5 ～ 8} \end{array}} \right. $

(5)

此时， $\tau {\text{ = }}\upsilon/c_s^2$ ，由流体黏性系数 $\upsilon$ 和声速 ${c_s} = 1/\sqrt 3$ 计算得到。

直接法求解的方程为式（1），而分步法则先采用θ方法将式（1）在时间区间[t,t + Δt]沿特征线方向积分，得到^[17]：

$ \begin{split} &{{f}_{\alpha }|}_{\left({\boldsymbol{x}}+{{\boldsymbol{e}}}_{\alpha }\Delta t,t+\Delta t\right)}-{{f}_{\alpha }|}_{({\boldsymbol{x}},t)} = {\frac{\theta }{{\tau }^{\prime }}\left({f}_{\alpha }^{\text{eq}}-{f}_{\alpha }\right)\Bigg|}_{({\boldsymbol{x}},t)}+\\ &\qquad{\frac{1-\theta }{{\tau }^{\prime }}\left({f}_{\alpha }^{\text{eq}}-{f}_{\alpha }\right)\Bigg|}_{\left({\boldsymbol{x}}+{{\boldsymbol{e}}}_{\alpha }\Delta t,t+\Delta t\right)} \end{split}$

(6)

上式 $ \theta $ 取为0.5以获得二阶时间精度， $ {\tau ^\prime }{\text{ = }}\tau /\Delta t $ ， $ \Delta t $ 为时间步长。为消除碰撞项的隐式，定义一个新的分布函数：

$  {\bar{f}}_{\alpha } = {f}_{\alpha }+\frac{{f}_{\alpha }-{f}_{\alpha }^{\text{eq}}}{2{\tau }^{\prime }} $

(7)

由上式可得：

$ \bar f_{\alpha} ^{{\text{eq}}} = f_{\alpha} ^{{\text{eq}}} $

(8)

因此，式（6）可以重写为：

$ {\bar f_{\alpha} }({\boldsymbol{x}},t) - {\bar f_{\alpha} }\left( {{\boldsymbol{x}} - {{\boldsymbol{e}}_{\alpha} }\Delta t,t - \Delta t} \right) = - {\left. {\frac{1}{{{\tau ^\prime } + 0.5}}\left( {{{\bar f}_{\alpha} } - \bar f_{\alpha} ^{{\text{eq}}}} \right)} \right|_{\left( {{\boldsymbol{x}} - {{\boldsymbol{e}}_{\alpha} }\Delta t,t - \Delta t} \right)}} $

(9)

对式（9）的求解分为两步，首先执行碰撞步^[17]：

$  {\hat{f}}_{\alpha } = {\bar{f}}_{\alpha }-\frac{{\bar{f}}_{\alpha }-{\bar{f}}_{\alpha }^{\text{eq}}}{{\tau }^{\prime }+0.5} $

(10)

接着求解以下的纯对流方程：

$ \frac{\partial {\hat{f}}_{\alpha }}{\partial t}+{{\boldsymbol{e}}}_{\alpha }\cdot\nabla {\hat{f}}_{\alpha } = 0 $

(11)

本文采用高精度FR格式对式（1）和式（11）进行空间离散。

1.2 空间离散

本文仅讨论四边形单元下FR格式，考虑以下的二维守恒律：

$  \frac{\partial w}{\partial t}+\nabla \cdot {\boldsymbol{F}}(w) = 0 $

(12)

其中， ${\boldsymbol{F}}{\text{ = }}({F_x},{F_y})$ 为通量矢量， $ {F_x} $ 、 $ {F_y} $ 分别为通量在x、y方向上的分量。

将空间划分为相互不重叠的四边形网格单元。为了在每个四边形单元采用统一的插值函数，简化计算过程，首先需要把不规则的单元映射到标准单元 $ [-1,1] \times [-1,1] $ 上，每个求解点在标准单元的坐标是一致的。物理空间全局坐标系采用（x，y），标准单元计算域局部坐标系采用（ $ \xi $ ， $\eta $ ），它们之间的可逆转换关系为：

$ \left( {\begin{array}{*{20}{l}} {x(\xi ,\eta )} \\ {y(\xi ,\eta )} \end{array}} \right) = \sum\limits_{i = 1}^K {{M_i}} (\xi ,\eta )\left( {\begin{array}{*{20}{l}} {{x_i}} \\ {{y_i}} \end{array}} \right) $

(13)

式中，K表示定义该网格单元的节点数量，M为坐标转换函数。根据坐标变换关系，守恒方程（12）在标准单元下的表达形式为：

$  \frac{\partial \tilde{w}}{\partial t}+\nabla \cdot\tilde{{\boldsymbol{F}}} = \frac{\partial \tilde{w}}{\partial t}+\frac{\partial {\tilde{F}}_{x}}{\partial \xi }+\frac{\partial {\tilde{F}}_{y}}{\partial \eta } = 0 $

(14)

式中，

$  \tilde{w} = \left|{\boldsymbol{J}}\right|w \text{，}\left[ {\begin{array}{*{20}{c}} {{{\tilde F}_x}} \\ {{{\tilde F}_y}} \end{array}} \right] = |{\boldsymbol{J}}|{{\boldsymbol{J}}^{ - 1}}\left[ {\begin{array}{*{20}{l}} {{F_x}} \\ {{F_y}} \end{array}} \right] $

(15)

其中雅克比矩阵 ${\boldsymbol{J}}$ 定义为：

$ {\boldsymbol{J}} = \frac{\partial (x,y)}{\partial (\xi ,\eta )} = \left(\begin{array}{c}{x}_{\xi },{x}_{\eta }\\ {y}_{\xi },{y}_{\eta }\end{array}\right) $

(16)

为了在单元内构造高阶的连续通量插值多项式，需要定义两类点集，即求解点和通量点。求解点本文应用Gauss-Legendre点，通量点为每个方向上的边界点。三阶FR方法的标准单元如图1所示。

以 $ \xi $ 方向的通量为例，通过两个一维拉格朗日插值基函数的张量积可以构造出二维的通量插值多项式：

$ {\tilde F_x}(\xi ,\eta ) = \sum\limits_{j = 1}^N {\sum\limits_{i = 1}^N {{{\tilde F}_{i,j}}} } \left[ {{l_i}(\xi ) \cdot {l_j}(\eta )} \right] $

(17)

其中，N为一维方向上的求解点个数， $ {\tilde F_{i,j}} $ 为单元内 $ \xi $ 方向第i个、 $ \eta $ 方向上第j个求解点处的通量，l为拉格朗日插值基函数。由拉格朗日插值基函数的性质可知，对于单元内每一个固定的 $ \eta $ ，沿着 $ \xi $ 方向问题可以简化为一维，即：

$ {\tilde F_x}(\xi ) = \sum\limits_{i = 1}^N {{{\tilde F}_i}} {l_i}(\xi ) $

(18)

以上的通量多项式并没有考虑单元间的作用，因此，它在单元边界是间断的。为了重新构造单元间连续的通量多项式，需要引入左右边界的公共通量 $ \left( {{{\tilde F}_x}} \right)_{i - \frac{1}{2}}^{{\text{com }}} $ 和 $ \left( {{{\tilde F}_x}} \right)_{i{\text{ + }}\frac{1}{2}}^{{\text{com }}} $ 进行修正，采用迎风数值通量Roe格式^[24]计算公共通量。定义修正后的通量多项式为 $ \tilde F_x^{{\text{c }}}\left( \xi \right) $ ，其表达式为：

$ \begin{split} \tilde F_x^{{\text{c }}}\left( \xi \right) =& {{\tilde F}_x}(\xi ) + \left[ {\left( {{{\tilde F}_x}} \right)_{i - \frac{1}{2}}^{{\text{com }}} - {{\tilde F}_x}( - 1)} \right]{g_{_{\rm{L}}}}(\xi )+ \\ & \left[ {\left( {{{\tilde F}_x}} \right)_{i{\text{ + }}\frac{1}{2}}^{{\text{com }}} - {{\tilde F}_x}(1)} \right]{g_{_{\rm{R}}}}(\xi ) \\ \end{split}$

(19)

其中 ${g_{_{\rm{L}}}}$ 和 ${g_{_{\rm{R}}}}$ 分别为左右边界的修正函数。这个函数是FR方法的关键，影响着格式的精度和稳定性。它是一个N阶多项式，即通量多项式的阶数等于一维方向求解点的数目。本文采用Vincent等提出的Vincent-Castonguay-Jameson-Huynh （VCHJ）格式修正函数，通过调整函数中的一个参数，使得FR方法可以恢复NDG、SD以及Huynh的 $ {g_2} $ 方案^[25]。并且，该格式被证明用来解线性对流方程以及对流扩散方程任意阶都稳定^[21,26]，具体形式参见文献[21]。为了表述方便，称采用N阶通量多项式的通量重构格子玻尔兹曼法为p^N FRLBM。

最终，对式（19）两边同时关于 $ \xi $ 求导，得到：

$ \begin{split} \frac{{{\rm{d}}\tilde F_x^{{\text{c }}}}}{{{\rm{d}}\xi }}\left( \xi \right) =& \frac{{{d}\tilde F_x^{\text{ }}}}{{{\rm{d}}\xi }}\left( \xi \right) + \left[ {\tilde F_{i - \frac{1}{2}}^{{\text{com }}} - {{\tilde F}_i}( - 1)} \right]\frac{{{\rm{d}}{g_{_{\rm{L}}}}}}{{{\rm{d}}\xi }}\left( \xi \right) +\\ &\left[ {\tilde F_{i + \frac{1}{2}}^{\text{com}} - {{\tilde F}_i}(1)} \right]\frac{{{\rm{d}}{g_{_{\rm{R}}}}}}{{{\rm{d}}\xi }}\left( \xi \right) \end{split} $

(20)

1.3 时间离散

如果采用显式时间离散格式求解方程（1），每个单元i的时间步长受限于CFL条件以及碰撞松弛时间，如下式所示：

$ \Delta {t_i} = \sigma \min \left( {\frac{{\sqrt {{V_i}} }}{{{\varLambda _i}}},2\tau } \right)\bigg /\left( {2N + 1} \right) $

(21)

其中 $ \sigma $ 、 $ {V_i} $ 和 ${\varLambda _i}$ 分别为CFL数、单元i的体积和对流谱半径。为摆脱松弛时间的限制，本文用两种方法，一种是Guo和Zhao提出的基于半隐格式的显式方法，具体实施方法如下。在时间区间 $ \left[ {{t_n},{t_{n + 1}}} \right] $ 对方程（1）进行积分可得：

$ f_{\alpha} ^{n + 1} - f_{\alpha} ^n + \Delta t{{\boldsymbol{e}}_{\alpha} }\cdot\nabla f_{\alpha} ^n = \Delta t\left[ {\theta \varOmega _{\alpha} ^{n + 1} + (1 - \theta )\varOmega _{\alpha} ^n} \right] $

(22)

其中， $\varOmega _{\alpha} ^{}$ 为右端碰撞项， $\theta $ 取为1/2。引入一个新的分布函数来消除碰撞项中的隐式：

$  {g}_{\alpha } = {f}_{\alpha }-\frac{1}{2}\Delta t{\varOmega }_{\alpha } $

(23)

式（22）重写为：

$ g_{\alpha} ^{n + 1} = f_{\alpha} ^n - \Delta t{{\boldsymbol{e}}_{\alpha} }\cdot\nabla f_{\alpha} ^n + \frac{1}{2}\Delta t\varOmega _{\alpha} ^n $

(24)

n+1时刻的宏观量可以对 $g_{\alpha} ^{n + 1}$ 求矩获得；接着， $f_{\alpha} ^{n + 1}$ 可以由下式求得：

$ f_{\alpha} ^{n + 1} = \frac{{g_{\alpha} ^{n + 1} + \dfrac{{\Delta t}}{{2\tau }}{{\left( {f_{\alpha} ^{{\text{eq}}}} \right)}^{n + 1}}}}{{1 + \dfrac{{\Delta t}}{{2\tau }}}} $

(25)

其中 ${\left( {f_{\alpha} ^{{\rm{eq}}}} \right)^{n + 1}}$ 可以由n+1时刻的宏观量代入式（2）求得。

另一种方法是隐式-显式Runge-Kutta法，其表达形式如下：

$ \begin{split} f_{\alpha} ^{(j)} =& f_{\alpha} ^n - \Delta t\sum\limits_{k = 1}^{j - 1} {{{\tilde m}_{jk}}} \left( {{{\boldsymbol{e}}_a}\cdot\nabla f_{\alpha} ^{(k)}} \right)- \hfill \\ & \Delta t\sum\limits_{k = 1}^j {{m_{jk}}} \frac{{f_{\alpha} ^{(k)} - f_{\alpha} ^{{\text{eq}}(k)}}}{\tau } \hfill \end{split} $

(26)

$ \begin{split} f_{\alpha} ^{n + 1} =& f_{\alpha} ^n - \Delta t\sum\limits_{j = 1}^r {{{\tilde n}_j}} \left( {{e_{\alpha} }\cdot\nabla f_{\alpha} ^{(k)}} \right) - \hfill \\ & \Delta t\sum\limits_{j = 1}^r {{n_j}} \frac{{f_{\alpha} ^{(k)} - f_{\alpha} ^{{\text{eq}}(k)}}}{\tau } \hfill \end{split} $

(27)

其中， $f_{\alpha} ^{(j)}$ 和 $ f_{\alpha} ^{{\text{eq}}(j)} $ 分别为中间第j步的分布函数和平衡态， $ {\tilde m_{jk}} $ 、 $ {m_{jk}} $ 、 $ {\tilde n_j} $ 、 $ {n_j} $ 为格式系数。式（26）等号右边第三项的隐式部分可以利用碰撞不变量的特性消除，格式的具体细节和不同精度的格式系数可以参见文献[16]。对方程（11）的时间推进则采用常用的三阶强稳定性Runge-Kutta法。最终的时间步长选取仅受CFL条件限制，如下式所示：

$ \Delta {t_i} = \sigma \min \left( {\frac{{\sqrt {{V_i}} }}{{{\varLambda _i}}}} \right)\bigg/\left( {2N + 1} \right) $

(28)

1.4 边界条件

本文数值算例中涉及的无滑移壁面及远场边界处理，采用非平衡态外推的方法^[27]，将边界上的分布函数写成平衡态和非平衡态之和的形式：

$ {\left( {{f_{\alpha} }} \right)_{\text{b}}} = {\left( {f_{\alpha} ^{{\text{eq}}}} \right)_{\text{b}}}{\text{ + }}{\left( {f_{\alpha} ^{{\text{neq}}}} \right)_{\text{b}}} $

(29)

其中，下标b表示边界，上标neq表示非平衡态。平衡态可以根据边界的宏观量直接求得，而非平衡态可以由求解点的非平衡态通过拉格朗日插值得到。一旦得到边界上的分布函数，就可以直接计算出边界通量，然后通过式（19）将边界的影响考虑进来。

2 数值算例 2.1 Taylor-Green涡

Taylor-Green涡问题被广泛用来研究非定常不可压流动算法的收敛性，这个问题的解析解为：

$ \left\{ {\begin{array}{*{20}{l}} {{u^*}(x,y,t) = - {U_0}\cos (2{\text{π}}x)\sin (2{\text{π}}y)\exp \left( { - 8{{\text{π}}^2}vt} \right)} \\ {{v^*}(x,y,t) = {U_0}\sin (2{\text{π}}x)\cos (2{\text{π}}y)\exp \left( { - 8{{\text{π}}^2}vt} \right)} \\ {{p^*}(x,y,t) = - \dfrac{{U_0^2}}{4}\big[\cos (4{\text{π}}x) + \cos (4{\text{π}}y)\big]\exp \left( { - 16{{\text{π}}^2}vt} \right)} \end{array}} \right. $

(30)

其中 $ {U_0} $ 为初始速度，为了减少格子玻尔兹曼的压缩性误差，取为0.001。计算域为一个边长为1的正方形，所有边界为周期性条件。雷诺数Re取为10，初始的速度场和压强场可以由式（30）取 $ t{\text{ = }}0 $ 获得，网格为均匀网格。x方向速度的相对误差L₂定义为：

$ \begin{split} {L_2}\left( t \right) =& {\left\| {u(t) - {u^*}(t)} \right\|_2} \\ = & \sqrt {\frac{{\displaystyle\sum\limits_{i = 1}^{{N_{{\text{total}}}}} {{{\left[ {{u_i}(t) - u_i^*(t)} \right]}^2}} }}{{\displaystyle\sum\limits_{n = 1}^{{N_{{\text{total}}}}} {{{\left[ {u_i^*(t)} \right]}^2}} }}} \end{split} $

(31)

其中，N_total为求解点总数， $ {u_i}(t) $ 为数值解， $ u_i^*(t) $ 为解析解，时间取 ${t_h} = \ln 2/\left( {8v{\text{π} ^2}} \right)$ 。为减少时间离散误差，时间步长取0.0002。本文比较了直接法与分步法的数值表现，其中直接法时间离散采用二阶IMEX。不同网格单元计算得到的误差和收敛精度阶见表1。从表中可知，误差会收敛到一个固定值，这是模型的压缩性误差导致的，这一现象在其他高精度OLBM中同样存在^[13-14]。在被压缩性误差影响前，直接法和分步法都能获得高阶收敛阶，并且直接法与分步法的误差及收敛阶几乎一致。以上结论是基于一个较小时间步长得到的，为了比较不同时间步长对两种方法结果的影响，图2给出了不同CFL数 $ \sigma $ 下、网格单元为16×16时，p³ FRLBM直接法和分步法的相对误差和计算时间，时间步长由式（28）算得。从图2可以看到，CFL数较小时两种方法误差几乎一致。但当CFL数增加，直接法的误差几乎不变，而分步法明显上升。这是由于分步法的对流方程可以看成是拉格朗日迁移的一种欧拉表达^[28]，标准LBM迁移是该对流方程不引入任何相位和幅值误差的精确解。当CFL数取0.7时，直接法的计算时间为19 s；而要取得与之接近的误差，分步法的CFL数最大取0.3，此时分步法计算时间为28 s。从这个结果可以看出，取得相似误差时，直接法能取更大的时间步长，效率更高。此外，直接法的最大CFL数可以取到0.85，而分步法最大只能取到0.7。因此，直接法的稳定性更好。

2.2 顶盖驱动方腔流

顶盖驱动方腔流因为其流动结构丰富、流场变化受雷诺数影响明显的特点，是不可压流算法的经典校核算例之一。在本算例中，计算域是一个边长L为1的方腔，顶盖驱动速度u₀为0.1，其他三个边为无滑移壁面，计算的雷诺数定义为Re = u₀L/ $\upsilon$ ，本文计算的四个稳态Re分别为3200、5000、7500和10000。计算收敛的判据为：

$ \frac{{\sqrt {\displaystyle\sum\limits_{{N_{{\rm{total}}}}} {\left[ {{{\left( {u_{}^{t + \Delta t} - u_{}^t} \right)}^2} + {{\left( {v_{}^{t + \Delta t} - v_{}^t} \right)}^2}} \right]} } }}{{\sqrt {\displaystyle\sum\limits_{{N_{{\rm{total}}}}} {\left[ {{{\left( {u_{}^{t + \Delta t}} \right)}^2} + {{\left( {v_{}^{t + \Delta t}} \right)}^2}} \right]} } }} < {10^{ - 8}} $

(32)

为了更好地捕捉方腔角落处的涡，采用非均匀网格：

$ x = \frac{1}{{2a}}[a + \tanh (\beta \xi )]\quad $

(33)

$ y = \frac{1}{{2a}}[a + \tanh (\beta \eta )] $

(34)

其中， $a = \tanh \beta$ ， $ \beta = 1.5 $ 。当Re = 3200，采用的网格单元数为30×30；当Re = 5000和7500，采用的网格单元数为40×40；当Re = 10000，采用的网格单元数为50×50。40×40非均匀网格单元的分布如图3所示。图4为p⁴ FRLBM计算得到的四个不同雷诺数下的流线。从图中可以看出，Re = 3200时右下角出现小涡，Re = 5000时左下角出现小涡。随着雷诺数的增加，底角三级涡越来越明显。由于FR格式的高阶精度和高分辨率特性以及非均匀网格的使用，这些流动细节被FRLBM清楚地捕获。图5定量比较了方腔中线的无量纲速度，其中左侧纵坐标和右侧纵坐标分别为x方向和y方向无量纲速度。本文结果和相关文献^[10,12,29]吻合良好，尽管本文采用了较少的网格单元。此外，本文还比较了4种时间离散方式在4个Re下收敛所能取到的最大CFL数，如图6所示。可以发现，基于半隐格式的显式方法和一阶IMEX所能取到的最大CFL数相同，为0.25。并且，随着雷诺数的增大，CFL数降低。而二阶IMEX和三阶IMEX即使在雷诺数较大时，CFL数依然可以取为1附近，且三阶IMEX最大CFL数略大于二阶IMEX。考虑到三阶IMEX要增加一步空间算子，二阶IMEX是较好的时间离散方式。

2.3 圆柱绕流

圆柱绕流是验证算法模拟复杂外形绕流的经典算例之一，本文计算了瞬态雷诺数Re = 200的情况。采用沿径向延展的非均匀O型贴体网格，网格单元为1次曲线单元，即式（13）中的K = 8。径向延展的变化由下式确定：

$ {r_j} = {r_1} + \left( {{r_1} - {r_0}} \right)\frac{{\tanh \left( {\beta {\eta _j}} \right)}}{{\tanh \beta }} $

(35)

其中，r₀为圆柱的半径，取为0.5；r₁为圆形计算域的半径，取为50； $ {\eta _j}{\text{ = }}j/{N_r} - 1 $ ，j = 0, 1, ··· , N_r；N_r 是径向网格单元的数量，取为60；周向单元数量取为40； $\; \beta $ = 3。计算域及网格单元的分布如图7所示。图8给出了瞬时流线和涡量等值线，从图中可以观察到卡门涡街的脱落。升力系数C_L和阻力系数C_D随时间的周期性变化见图9，由图可以观察到C_D的周期约为C_L的一半。此外，将本文计算得到的C_D、C_{L ，}以及斯特劳哈尔数St和相关文献^[7,11,30,31]做了定量比较，如表2所示，可以发现吻合良好。斯特劳哈尔数的定义为 $S t = D/{u_0}{T_p}$ ，其中 $ D $ 为圆柱直径， $ {u_0} $ 为来流速度， $ {T_p} $ 为C_L随时间变化的峰值周期。因此，本文的方法模拟复杂外形绕流是可靠的。

3 结　论

OLBM可以改善标准LBM的几何适应性，本文回顾了不可压流动高精度OLBM的发展过程。相比基于N-S方程的高精度不可压流算法，基于DVBE的高精度算法要容易实现得多。本文基于高精度通量重构格式发展了一种通量重构格子玻尔兹曼方法，并在此框架下比较了直接求解DVBE和先碰撞、再求解纯对流方程这两种方法的数值表现。研究发现，在小时间步长下，这两种方法误差接近，采用高精度格式都能获得高阶的收敛阶。然而，当时间步长增大时，分步法的误差出现明显上升，而直接法误差几乎不变。数值测试显示，在取得相似误差时直接法能取更大的时间步长，运算周期较短、计算效率更高。此外，采用二阶隐式-显式Runge-kutta推进的直接法，稳定性优于分步法。接着，通过顶盖驱动方腔流算例，比较了4种碰撞项隐式处理的时间推进方法在不同雷诺数下所能取到的最大CFL数。数值结果表明，二阶IMEX是较好的时间离散方法。然而，尽管IMEX能使时间步长的选取摆脱松弛时间的限制，但由于对流项依然显式处理，在模拟壁湍流问题时较小的边界层网格厚度同样会严重限制时间步长。因此，发展低存储、全隐式的时间推进方法是OLBM的重要发展方向之一，目前已有学者开展了一些工作^[32-33]。本文采用的碰撞模型是单松弛，为了增强稳定性，也可以采用更先进的碰撞模型，比如多松弛等。此外，高精度OLBM还可以直接拓展到湍流、多相流、热流等应用中^[34-35]，目前相关工作正在开展中。