摘要:在传统中, 格子Boltzmann离散模型的数值表现影响因素分析主要关注于平衡分布的矩精度, 而模型的正值性(作为粒子分布描述, 分布函数需恒为正)则仅作为模型的附属属性, 用于计算时工况约束. 随着部分高斯-厄米特求积公式离散方案的提出, 模型正值性被发现是独立于矩精度的模型属性, 可以通过格子速度调整. 研究人员推测平衡分布正值性对格子Boltzmann方法的数值表现存在显著影响, 可以通过改善平衡分布正值性改善模型数值表现. 相比提升模型矩精度方案, 正值性改善方案具有计算量优势. 然而鉴于高阶模型边界处理的缺失, 相关推测并未得到具体数值计算证实. 本文采用周期边界的Taylor-Green涡算例, 回避了边界处理问题, 详细分析了正值性对数值表现的影响, 包括平衡分布正值区域内计算精度稳定性以及模型平衡分布正值区域大小对计算影响, 并与矩精度影响进行对比. 计算结果显示, 模型正值区域内计算精度并不恒定, 随着工况靠近正值区域上界, 计算精度下降, 但总体上均具备较好的精算精度. 模型数值表现同时受到矩精度与模型正值性影响, 矩精度影响主要体现在模型是否满足Galilean不变性上, 对于满足Galilean不变性模型, 其数值表现则取决于模型正值性. 基于此, 本文认为通过改善模型正值性提升格子Boltzmann方法数值表现是切实可行的方案, 并推荐基于满足Galilean不变性条件下选择具有最宽正值区域的模型, 而不必执着于模型矩精度. 另外从本文的数值结果来看, 高阶模型模型的数值表现均好于经典D2Q9模型, 特别是D2H3-2模型, 是文中涉及模型的最优者, 值得进一步深入研究. 总体而言, 通过数值分析首次系统性梳理了模型正值性对计算影响, 并与矩精度进行对比分析. 本文证实了正值性对计算的影响, 为离散模型选择和改进提供了新的方向.
关键词: 格子Boltzmann/
离散模型/
矩精度/
正值性

English Abstract

全文HTML

--> --> -->

最近几十年, 格子Boltzmann方法在计算流体领域受到了广泛的关注. 作为一种介观尺度计算方案, 格子Boltzmann方法既能提供相较于宏观计算方法更为深入的视角, 又能避免纯粒子尺度计算方法所需要的大量计算资源, 因此在诸多传统方法难以适用的复杂流体领域得到了广泛应用, 如多组分多相流^[1-4], 微纳米尺度流^[5-7]、多孔介质流^[6,8]、粒子悬浮流^[9,10]、血液流^[11,12]、湍流^[13-15]等. 格子Boltzmann方法描述的是一群虚构粒子在流场内迁移和碰撞的过程. 流场内粒子信息通过粒子分布函数来描述, 而流场的宏观信息通过对粒子分布函数做相应的速度矩积分得到. 格子Boltzmann方法的控制方程可以通过对BGK近似下的Boltzmann方程沿特征线积分得到^[16-18], 而具体的算法实现则还需对得到的特征线积分在速度空间上进行离散^[17,18], 如常见的速度离散模型D2Q9等. 然而低阶矩精度速度离散模型在还原宏观控制方程式的过程中会引入误差, 如D2Q9模型还原动量守恒方程时会引入额外黏性项, 导致方程不满足Galilean不变性^[19-21]. 为解决该问题, 研究者们提出了高阶速度离散模型, 如Shan^[22,23]和Shim^[24]的工作. 然而由于方案实现的复杂性, 高阶速度模型数量非常有限, 导致对于高阶模型数值表现的理解集中于模型的矩精度.
最近Ye等^[25]提出了部分高斯-厄米特求积公式(partial Gauss-Hermite quadrature, pGHQ), 显著简化了离散模型的构建. Ye等^[25]用该方案系统性地搜索了矩精度介于{3—7}离散速度在$ [-10,10] $范围内的可行离散模型, 结果显示可行模型的数量非常丰富, 远超预期, 尤其是存在大量精度相同但是离散速度不同的模型. 这些模型从宏观角度分析是等价的, 即用Chapman-Enskog展开可以还原出相同的宏观方程. 但是在微观上, 模型有着明显的区别, 特别是平衡分布的正值区域(所有分布为正值的宏观速度区域). 为方便讨论, 本文用正值性指代分布正值相关性质包括正值区域等, 且在后续讨论中不区分两者使用. 相对于传统理解中作为模型用于计算时工况的约束条件, 在pGHQ方案下, 平衡分布正值性可以看作是独立于矩精度的模型属性, 可以通过离散速度调整. 鉴于负值平衡分布违背格子Boltzmann方法的物理原则, 即粒子分布不存在负值, Ye等^[25]认为离散模型的正值性将对模型的数值表现有直接影响, 可以通过选择模型离散速度来改善计算. 对比通过提升模型矩精度改善格子Boltzmann方法数值表现的传统构想, 调整模型离散速度设想在计算量上存在优势. 然而鉴于高阶模型的具体实现需要设计对应的边界处理, Ye等^[25]文章中并没给出具体数值验证. 本文选择可以采用周期性边界处理的Taylor-Green涡算例, 回避了高阶离散模型边界处理匮乏的问题. 通过详细的数值实验, 文章分析了平衡分布正值性对离散模型数值稳定性和精度的影响, 并与模型矩精度影响进行对比, 验证通过提升模型正值性来改善格子Boltzmann方法数值表现的可行性.
本文首先回顾Ye等^[25]提出的pGHQ离散模型构造方案, 并选择了部分离散模型用作具体计算. 然后本文用选择的模型计算Taylor-Green涡, 分析平衡分布正值性对格子Boltzmann方法计算稳定性和精度的影响, 并与模型矩精度影响进行对比, 验证通过提升模型正值性来改善格子Boltzmann方法数值表现的可行性. 最后总结本文工作.

pGHQ速度离散方法基于Hermite展开理论^[26], 即平衡分布可以通过对Maxwell分布以Hermite多项式形式展开构造. 以一维为例, 对于m阶精度的平衡分布, 其分布函数可以通过如下方式构造:

$ f_\alpha ^{\rm eq} = {{{w}}_\alpha }\rho \sum\limits_{i = 0}^m {{{ H}_i}\left( {\frac{{{v_\alpha }}}{{\sqrt 2 {c_{\rm s}}}}} \right)\frac{{{u^i}}}{{{{\left( {\sqrt 2 {c_{\rm s}}} \right)}^i}i!}}} , $

式中$ \rho $和u分别为密度和宏观流速, $ {v_\alpha } $和$ {w}_\alpha $为离散格子速度和对应的权重系数, $ 1/\sqrt 2 {c_{\rm s}} $在Hermite展开理论中被称为无量纲网格常数c, 而$ {{ H}_i}\left( x \right) $则为i阶物理Hermite多项式, 表1罗列了前5阶物理Hermite多项式的具体表达. 无量纲网格常数c, 离散格子速度$ {v_\alpha } $和对应的权重系数$ {w}_\alpha $是通过求解对应的$ 0-2m $阶数值积分方程得到,

Hermite 多项式阶数 The degree of Hermite polynomial	多项式表达式 Hermite polynomial
0	${ { H}_0}\left( x \right) = 1$
1	${ { H}_1}\left( x \right) = 2x$
2	${ { H}_2}\left( x \right) = 4x^2-2$
3	${ { H}_3}\left( x \right) = 8x^3-12x$
4	${ {H}_4}\left( x \right) = 16x^4-48x^2+12$

表1前五阶物理Hermite多项式
Table1.The first five physicists′ Hermite polynomials.

$ {I^k} = \int {\frac{1}{{\sqrt {\text{π}} }}{{\rm e}^{ - {\xi ^2}}}{\xi ^k} {\rm d}\xi } = \sum\limits_{\alpha = 0}^{q - 1} {{{{w}}_\alpha }\xi _\alpha ^k},\quad k \in [0,2m], $

式中$ \xi_\alpha = v_\alpha c $. 方程(2)是高阶非线性方程组, 为减少其求解难度, Shan^[22,23]和Shim^[24]预先定义了对称格子速度集$ \left\{ {0, \pm {v_1},\cdots} \right\} $, 使得方程组转变为仅包含无量纲网格常数c和权重系数$ \{{w}_\alpha\} $未知量的方程, 并假定m阶矩精度平衡分布构造需要$ 2m-1 $个格子速度. 至此离散速度模型构造转变为构造$ \left\{ {0, \pm {v_1},\cdots, \pm {v_{m - 1}}} \right\} $并求解其对应方程组(2). 尽管Shan和Shim等通过预定义格子速度集$ \left\{ {0, \pm {v_1},\cdots, \pm {v_{m - 1}}} \right\} $显著简化了离散速度模型构造, 但是求解代入格子速度集后的方程组(2)依旧繁琐而复杂, 无法大规模验证可行的格子速度集. 因此现有的高阶可行离散模型数量十分有限, 这阻碍了对于格子Boltzmann速度离散模型的进一步研究. 另外对于模型矩精度与格子速度数量之间的关系也没经过系统而深入地分析. 因此在Shan和Shim的研究中, 其重点更偏重于突破高阶离散模型的缺失, 并未形成高阶模型系统性认识. 这使得对于高阶模型的认知依旧局限于传统“模型阶数决定模型数值表现”的猜测, 即模型阶数越高数值计算表现越好.
Ye等^[25]通过改造Gauss-Hermite求积公式, 给出了方程(2)的等价求积公式, 即pGHQ求积公式. 根据pGHQ求积公式, 互不相同的q节点$ \xi_\alpha $和其权重系数$ {w}_\alpha $若要满足方程组(2), 只需保证其节点多项式$ {W_q}\left( \xi \right) $基于Hermite多项式表述时, 仅涉及$ 2m-q+1 $到q阶之间的Hermite多项式,

$ {W_q}\left( \xi \right)\!=\! \prod\limits_{\alpha = 0}^{q - 1} {\left( {\xi - {\xi _\alpha }} \right)}\!=\! \sum\limits_{i = 2m - q + 1}^{2m - 1} {{A_i}{{ H}_i}\left( \xi \right)} + \frac{1}{{{2^q}}}{{ H}_q}\left( \xi \right) , $

而对应的权重系数计算公式则为

$ {{{w}}_\alpha } = \frac{1}{{\sqrt {\text{π}} }}\int {{{\rm e}^{ - {\xi ^2}}}\prod\limits_{{\beta = 0},\; {\beta \ne \alpha }}^{q - 1} {\frac{{\xi - {\xi _\beta }}}{{{\xi _\alpha } - {\xi _\beta }}}} {\rm d}\xi }. $

完整的推导可参考文献[25]. 不同于(2)式同时包含节点$ \xi_\alpha $和权重系数$ {w}_\alpha $, (3)式仅为节点$ \xi_\alpha $的函数. 因此当给定格子速度集$ \{v_\alpha\} $时, 可以根据(3)式直接得到网格常数c的约束方程, 其具体流程是将格子速度集代入其节点多项式并展开多项式,

$ {W_q}\left( \xi \right) = \prod\limits_{\alpha = 0}^{q - 1} {\left( {\xi - {v_\alpha }c} \right)} = \sum\limits_{i = 0}^q {{b_i}{\xi ^i}}, $

其中$ b_i $是网格常数c的多项式, 引入次方项$ \xi^i $与Hermite多项式的关系式,

$ {\xi ^i} = \frac{{i!}}{{{2^i}}}\sum\limits_{l = 0}^{\left\lfloor {i/2} \right\rfloor } {\frac{1}{{l!\left( {i - 2l} \right)!}}{{ H}_{i - 2l}}\left( \xi \right)}, $

则(5)式可以表达为

$ {W_q}\left( \xi \right) = \sum\limits_{i = 0}^q {{b_i}{\xi ^i}} = \sum\limits_{i = 0}^{q - 1} {{a_i}{{ H}_i}\left( \xi \right)} + \frac{1}{{{2^q}}}{{ H}_q}\left( \xi \right). $

而根据pGHQ (3)式, 要使得格子速度满足m阶矩精度离散模型构造, 其节点多项式基于Hermite多项式表述时, 应当仅涉及$ 2m-q+1 $到q阶之间的Hermite多项式, 即$ 0 $到$ 2m-q $阶Hermite多项式系数为$ 0 $,

$ \left\{ {{a_i} = 0|i \leqslant 2m - q} \right\}. $

而系数$ a_i $和$ b_i $一样, 是仅包含网格常数c的多项式. 因此(8)式事实上是给定格子速度集$ \{v_\alpha\} $的网格常数c约束方程组. 求解(8)式, 若方程组有解, 则将求解的网格常数以及格子速度代入(4)式, 即可得到对应的权重系数. 将网格常数, 格子速度以及权重系数代入(1)式, 即可得到m阶矩精度离散模型.
在Hermite展开理论中, 高维模型可以通过对一维模型做张量乘积得到. 以二维为例, 其平衡函数构造即为

$ \begin{split}f_\alpha ^{{\rm{eq}}} =\; & \rho {w_{\alpha x}}{w_{\alpha y}}\left( {\sum\limits_{i = 0}^m {{{{H}}_i}\left( {\frac{{{v_{\alpha x}}}}{{\sqrt 2 {c_{\rm{s}}}}}} \right)\frac{{u_x^i}}{{{{\left( {\sqrt 2 {c_{\rm{s}}}} \right)}^i}i!}}} } \right)\\ & \times \left( {\sum\limits_{i = 0}^m {{{{H}}_i}\left( {\frac{{{v_{\alpha y}}}}{{\sqrt 2 {c_{\rm{s}}}}}} \right)\frac{{u_y^i}}{{{{\left( {\sqrt 2 {c_{\rm{s}}}} \right)}^i}i!}}} } \right),\end{split}$

对应的格子速度为

$ \left( {{v_{\alpha x}},{v_{\alpha y}}} \right) = \left\{ {\left( {{v_\beta },{v_\gamma }} \right)|{v_\beta },{v_\gamma } \in \left\{ {{v_\alpha }} \right\}} \right\}. $

Shan^[22,23]和Shim^[24]的离散方案, pGHQ速度离散方法显著简化了离散模型构建过程. 其算法的简洁和通用性让系统性研究格子速度集成为可能. Ye等^[25]用该方案系统性地搜索了矩精度分别为{3—7}，离散速度在$ [-10,10] $范围内的可行离散模型, 得到了大量等矩精度的离散模型. 这些模型矩精度相同, 但是正值性存在显著差别. 鉴于负值平衡分布违背格子Boltzmann方法的物理原则, 即粒子分布不存在负值, Ye等^[25]认为模型正值性会影响模型的数值表现, 并推测可以通过选择格子速度集改善模型正值性的方式改善格子Boltzmann方法的计算稳定性和精度.
为验证该推论的可行性, 本文选择了表2所列的离散模型, 模型以维度n和速度矩精度m来表示, 即DnHm. 文章主要从如下几个方面验证正值性影响: 1)平衡分布正值性对于数值表现影响, 是否在正值区域内, 模型都能保持计算的稳定性和相应的精度; 2) 平衡分布正值性与模型矩精度对数值表现影响对比; 3)通过调整平衡分布正值性改善数值表现可行性. 鉴于非对称网格对于平衡分布的正值性影响主要在于平移正值区域, 其在理论上更适合具有主流速度的计算问题. 而本文采用的Taylor-Green涡流场是对称性, 其主流速度为0, 因此并不适合针对于非对称网格的分析. 而从本文研究目标来看, 单纯采用对称离散模型足够确认相关结果. 因此本文在模型选择上并不涉及非对称离散模型. 同时为验证Hermite展开理论的实用性, 本文还额外采用了经典D2Q9模型作为比较,

Model name	Discrete velocityset {$v_\alpha$}	Lattice constant c	Weights {${w}_\alpha$}
D2H2	{$0, \pm 1$}	$1.2247\times10^0$	{$6.6667\times10^{-1}$, 1.6667$\times10^{-1}$}
D2H3-1	{$0, \pm 1, \pm 3$}	$5.5343\times10^{-1}$	{$7.4464\times10^{-2}$, $4.1859\times10^{-1}$, $4.4182\times10^{-2}$}
D2H3-2	{$0, \pm 2, \pm 5$}	$3.4420\times10^{-1}$	{$3.1044\times10^{-1}$, $3.0997\times10^{-1}$, $3.4812\times10^{-2}$}
D2H4	{$0, \pm 1, \pm 2, \pm 3$}	$8.4639\times10^{-1}$	{$4.7667\times10^{-1}$, $2.3391\times10^{-1}$, $2.6938\times10^{-2}$, $8.1213\times10^{-4}$}
D2H5	{$0, \pm 1, \pm 2, \pm 3, \pm 5$}	$4.7940\times10^{-1}$	{$1.6724\times10^{-1}$, $3.0315\times10^{-1}$, $5.3303\times10^{-2}$, $5.7922\times10^{-2}$, $2.0013\times10^{-3}$}
D2H6	{$0, \pm 1, \pm 2, \pm 3, \pm 4, \pm 5$}	$6.8590\times10^{-1}$	{$3.8694\times10^{-1}$, $2.4178\times10^{-1}$, $5.8922\times10^{-2}$, $5.6153\times10^{-3}$, $2.0652\times10^{-4}$, $3.2745\times10^{-6}$}

表2格子Boltzmann离散模型. 为方便展示, 表格罗列的是用于对应构造高阶模型的一维模型参数. 实际计算中需根据(9)式和(10)式对表格中的参数进行张量构造. 另外需注意的是表格中罗列的是网格常数c, 其与格子声速$ c_{\rm s} $存在如下换算关系$ c = 1/\sqrt{2}c_{\rm s} $
Table2.Discrete model description of lattice Boltzmann method. For the seeking of space saving, the parameters illustrated in this table are of the corresponding unidimensional models. In numerical implementations, they require tensor product illuminated in Eq. (9) and Eq. (10). It should be noted that the table lists the lattice constant c instead of the lattice sonic speed $ c_{\rm s} $, which can be expressed as $ c = 1/\sqrt{2}c_{\rm s} $.

$ f_\alpha ^{\rm eq} = {{{w}}_\alpha }\rho \left( {1 + \frac{{{{ v}_\alpha } \cdot u}}{{c_{\rm s}^2}} + \frac{{{{\left( {{{ v}_\alpha } \cdot u} \right)}^2}}}{{2c_{\rm s}^4}} - \frac{{{u^2}}}{{2c_{\rm s}^2}}} \right), $

其中格子速度, 权重系数以及格子声速$ c_{\rm s} $的取值与D2H2模型相同.
对于DnHm类型的模型, 其正值区域是比较容易分析的, 即只需分析对应一维模型的所有分布保持正值性的流体速度区域, 高维模型的正值区域即为该区域的张量乘积结果. 对于本文涉及的$ \{0,\pm v_1,\cdots, $$\pm v_{m-1}\}$类型一维模型, 其正值区域具有对称性. 经过张量乘积后, 该对称性在高维下依旧得到保留, 因此本文以一维正值区域上限来指代模型正值区域, 符号记为$ U_{\rm pos} $. 以本文D2H2模型对应的一维模型为例, 其平衡分布在$ [-0.82,0.82] $之间均为正值, 超过该范围后则分布出现负值. 而D2H2模型是一维模型的张量乘积结果, 因此流体速度在$ [-0.82,0.82]\times[-0.82,0.82] $区间内, D2H2模型均为正值, 即D2H2模型正值区域为$ [-0.82,0.82]\times $[–0.82, 0.82], 其$ U_{\rm pos} = 0.82 $. 通过同样分析可以得到本文其他DnHm型的正值区域. 比较特殊的是D2Q9模型, 不同于DnHm类型模型的正值区域为一维张量乘积结果, D2Q9模型下不同方向的流体速度是相互关联的. 为方便分析, 取D2Q9模型正值区域中最大$ U_{\rm pos} $, 即在$ [-U_{\rm pos},U_{\rm pos}]\times[-U_{\rm pos},U_{\rm pos}] $上满足D2Q9正值性的最大$ U_{\rm pos} $值. 本文涉及模型的正值区域如图1所示.
图 1 格子Boltzmann离散模型正值性分析　虚线框表示模型正值区域, 其涵盖区域为$ [-U_{\rm pos},U_{\rm pos}]\times[-U_{\rm pos},U_{\rm pos}] $, 模型$ U_{\rm pos} $标注在图中模型名称之后
Figure1. Positivity analysis of lattice Boltzmann discrete models. The dash squares denote the positive areas of each lattice Boltzmann discrete model, which are constructed as $ [-U_{\rm pos},U_{\rm pos}]\times[-U_{\rm pos},U_{\rm pos}] $. The $ U_{\rm pos} $ values of each model are annotated after the model names.

两维Taylor-Green涡作为等温不可压Navier-Stokes方程组的有限几组精确解, 是流体计算中常用的经典算例, 其控制方程为:

$ \begin{split}& \frac{{\partial u}}{{\partial x}} + \frac{{\partial v}}{{\partial y}} = 0,\\& \frac{{\partial u}}{{\partial t}} + u\frac{{\partial u}}{{\partial x}} + v\frac{{\partial u}}{{\partial y}} = - \frac{1}{\rho }\frac{{\partial p}}{{\partial x}} + \nu \left( {\frac{{{\partial ^2}u}}{{\partial {x^2}}} + \frac{{{\partial ^2}u}}{{\partial {y^2}}}} \right),\\& \frac{{\partial v}}{{\partial t}} + u\frac{{\partial v}}{{\partial x}} + v\frac{{\partial v}}{{\partial y}} = - \frac{1}{\rho }\frac{{\partial p}}{{\partial y}} + \nu \left( {\frac{{{\partial ^2}v}}{{\partial {x^2}}} + \frac{{{\partial ^2}v}}{{\partial {y^2}}}} \right).\end{split}$

其解为

$ \begin{split} & u = { - {u_0}\sqrt {{k_y}/{k_x}} \cos \left( {{k_x}x} \right)\sin \left( {{k_y}y} \right){{\rm e}^{ - t/{T_{\rm d}}}}},\\ & v = {{u_0}\sqrt {{k_x}/{k_y}} \sin \left( {{k_x}x} \right)\cos \left( {{k_y}y} \right){{\rm e}^{ - t/{T_{\rm d}}}}},\\ & p = {{p_0} \!-\! {\rho _0}\frac{{u_0^2}}{4}\!\bigg(\!{\frac{{{k_y}}}{{{k_x}}}\cos \left( {2{k_x}x} \right) \!+\! \frac{{{k_x}}}{{{k_y}}}\cos \left( {2{k_y}y} \right)} \!\!\bigg){{\rm e}^{ - 2t/{T_{\rm d}}}}}, \end{split} $

式中$ \rho_0 $为流体常物性密度; $ p_0 $为参考压力; $ u_0 $为Taylor-Green涡特征速度; k_x, k_y为周期因子; T_d为衰减特征时间，其计算式为$ 1/{T_{\rm d}} =$ $ \nu \left( {k_x^2 + k_y^2} \right) $. 为方便讨论, 本文将周期因子恒定选择为$ k_x = k_y = 1 $. 根据两维Taylor-Green涡的解, 当流场初始速度和压力分布满足:

$ \begin{split} & u ={ - {u_0}\cos x \sin y},\\ &v = {{u_0}\sin x \cos y},\\ &p = {{p_0} - {\rho _0}\frac{{u_0^2}}{4}\left[ {\cos \left( {2x} \right) + \cos \left( {2y} \right)} \right]}. \end{split} $

其速度场和压力场分别以指数$ \exp(-t/T_{\rm d}) $和$ \exp(-2t/T_{\rm d}) $衰减, 且速度场和压力场在流场内分别以$ 2{\text{π}} $和$ {\text{π}} $为周期成周期分布. Taylor-Green涡解十分有标识性, 统一的指数衰减形式非常适合用于测试算法的稳定性. 而周期分布避免了边界处理, 适用于验证目前边界处理缺乏的高阶离散模型.
本文用格子Boltzmann计算了流域为$ [0,2{\text{π}}]\times $$[0,2{\text{π}}] $的Taylor-Green涡演变. 流体的运动黏性$ \nu $为$ 0.1\;\rm{m^2/s} $, 特征速度$ u_0 $为$ 1 \; \rm{m/s} $. 网格密度为$ 200\times 200 $. 模拟采用单松弛因子格子Boltzmann方程,

$\begin{split} & {f_\alpha }\left( {{{r}} + \Delta t{{ v}_\alpha },t + \Delta t} \right) \\ =\; & {f_\alpha }\left( { r},{{t}} \right) - \frac{1}{\tau }\left[ {{f_\alpha }\left( { r},{{t}} \right) - f_\alpha ^{\rm eq}\left( { r},{{t}} \right)} \right],\end{split}$

式中$ {f_\alpha }\left( { r},{{t}} \right) $, $ {f_\alpha ^{\rm eq}\left( { r},{{t}} \right)} $分别表示t时刻$ { r} $节点上离散速度$ { v}_\alpha $对应的粒子分布和平衡分布, $ {f_\alpha }({ r} + \Delta t{{ v}_\alpha }, $$t + \Delta t) $表示经过一个计算步, $ t + \Delta t $时刻$ {{ r} + \Delta t{{ v}_\alpha }} $节点上离散速度$ { v}_\alpha $对应的粒子分布, $ \tau $则是松弛因子, 其计算式为

$ \tau =\frac{{\rm{1}}}{{\rm{2}}}+\frac{{\nu \Delta t}}{{c_{\rm s}^2\Delta {L^2}}} , $

其中$ \Delta t $和$ \Delta L $分别为计算时间步长和网格长度, $ \nu $为流体运动黏性, $ c_{\rm s} $为格子声速. 流场的初始分布为:

$ \begin{split} & {u_{L{\rm B}}} = - {u_{L{\rm B} ,0}}\cos x\sin y ,\\ & {v_{L{\rm B}}} = {u_{L{\rm B} ,0}}\sin x \cos y, \end{split} $

其中$ u_{L{\rm B} ,0} $为$ u_0 $在格子Boltzmann体系下的值, $ u_{L{\rm B} ,0} = u_0\Delta t/\Delta L $. 不同于初始流场引入, 初始压力场的引入需经过转换处理. 本文根据压力和密度关系式$ p = \rho c_{\rm s}^2 $, 通过设置初始密度场来引入压力初始条件,

$ \rho ={\rho _{{p_0}}} - {\rho _0}\frac{{u_{L{\rm B} ,0}^2}}{{4c_{\rm s}^2}}\left( {\cos \left( {2x} \right) + \cos \left( {2y} \right)} \right) . $

为避免初始密度出现负值, 本文统一取$ {\rho _{{p_0}}} = 10 $, $ {\rho _{0}} = 1 $. 从Taylor-Green涡的设置可以看出, 要验证平衡分布正值性影响, 可以通过调整$ u_{L{\rm B} ,0} $来实现, 即通过调整计算步长来改变$ u_0 $在格子Boltzmann体系下的值, 人为控制流场初始时刻平衡分布的正值性. 根据模型正值区域分析结果(图1), 本文设计了如表3中所列的计算工况.

Case name	Step time/s	$u_{L{\rm B} ,0}$
a	$1.570796\times10^{-2}$	0.5
b	$3.141593\times10^{-2}$	1.0
c	$4.712389\times10^{-2}$	1.5
d	$6.283185\times10^{-2}$	2.0

表3算例设计
Table3.Case design for lattice Boltzmann method

根据Taylor-Green涡的初始条件设置, 随着$ u_{L{\rm B} ,0} $增大, 格子Boltzmann的初始分布将开始出现负值, 而具体的负值情况则取决于平衡分布模型. 为方便分析和展现具体分布正值情况, 本文设计了如下正值指标:

$ {\kappa} = \frac{{\displaystyle\sum {\left| f \right|} }}{{\displaystyle\sum f }}=\frac{{\displaystyle\sum {\left| f \right|} }}{\rho }. $

当分布中存在负值时, $ \kappa $值将大于1, 而不存在负值时则$ \kappa $值恒为1. $ \kappa $偏离1的幅值可以看作分布正值性破坏的严重程度. 该指标通过对$ \rho $归一, 消除了Taylor-Green涡流场内不同密度对正值性的影响. 图2给出了设计计算工况下, 各模型初始分布的$ \kappa $值情况. 图中额外标注了对应工况下模型正值区域值$ u_{\rm pos} $与设计$ u_{L{\rm B} ,0} $的差值, 记作$ \Delta u $,
图 2 各格子Boltzmann模型在Taylor-Green涡设计计算工况下初始时刻计算流域$ [0,2{\text{π}}]\times[0,2{\text{π}}] $内$ \kappa $值色阶图　图中列对应格子Boltzmann模型, 行则对应表3中计算工况; 每个色阶图上方的数值$ \Delta u $值, 即正值区域值$ u_{\rm pos} $设计计算工况$ u_{L{\rm B} ,0} $差值;为方便观察正值性随计算工况遭到破坏的变化过程, 图中对$ \kappa $值做了截断处理, 仅区别$ [1.0,2.0] $范围变化
Figure2. The initial $ \kappa $ filled contours of designed Taylor-Green vortex simulations for each lattice Boltzmann model in the fluid domain $ [0,2{\text{π}}]\times[0,2{\text{π}}] $. Each column is a simulation set of a lattice Boltzmann model with different designed cases. Each row is a set of a simulation case with different lattice Boltzmann models. The values in the contour panels are the values of $ \Delta u $, i.e. the difference between the $ u_{\rm pos} $ and the designed $ u_{L{\rm B} ,0} $. For the sake of identifying the positivity violation process of a lattice Boltzmann model, the filled contours truncate the range of $ \kappa $, limiting into $ [1.0,2.0] $.

$ \Delta u = u_{\rm pos} - u_{L{\rm B} ,0}. $

用以辅助理解设计计算工况对模型正值区域的偏离程度. $ \Delta u $值为正, 表示计算工况完全位于模型正值区域内, $ \Delta u $值为负则表示计算工况超过正值区域, 部分网格出现负值分布. 另外需注意的是, 图2对$ \kappa $值的范围做了截断处理, 即仅区别$ [1.0, 2.0] $范围内变化, 以方便观察分布正值性随设计工况而遭到破坏的变化过程. $ \kappa $值图显示, 当$ \Delta u $值为负时, 分布的正值性遭到破坏, $ \kappa $值大于1, 如D2H2模型在b, c, d工况情况下所示. 但是对于给定负值$ \Delta u $, 平衡分布正值性遭到的破坏程度取决于具体模型. 以D2H2和D2H6模型为例, D2H2在$ \Delta u $值为$ -0.18 $时, 初始$ \kappa $分布出现了明显偏离$ 1 $值; 而D2H6模型在$ \Delta u $值为$ -1.02 $时, 初始$ \kappa $分布偏离$ 1 $值的幅度依旧小于前述情况.
对于模型数值表现分析, 本文主要考量计算的流场分布以及误差演化. 流场分布分析通过比较Taylor-Green涡初始流场强度衰减一半时(t = $ 3.4657359\; \rm{s} $)模型计算结果与理论解的速度平方和$ U^2 $等高线图实现,

$ U^2 = u^2+v^2. $

鉴于速度平方和$ U^2 $在流域内成周期性分布, 本文仅截取周期区域$ [0.3775,0.6225] {\text{π}} \times[0.3775,0.6225] {\text{π}} $内图像(见图3). 误差演化分析则关注模型计算结果与理论解误差在$ t = 3.4657359 \;\rm{s} $时间内随时间的演变(图4), 其计算公式为
图 3 Taylor-Green涡半衰时刻$ t = 3.4657359\; \rm{s} $, $ [0.3775,0.6225]{\text{π}} \times[0.3775,0.6225]{\text{π}} $流域内, 不同计算工况下理论解(Theoretical) 与各格子Boltzmann模型的速度平方和$ U^2 $等高线图; 图中列对应格子Boltzmann模型, 行则对应表3中计算工况
Figure3. The $ U^2 $ contour lines of Taylor-Green vortex in fluid domain $ [0.3775,0.6225]{\text{π}} \times[0.3775,0.6225]{\text{π}} $ at half-value decay time $ t = 3.4657359 \; \rm{s} $, including the theoretical solution (denoted as “Theoretical”) and the numerical results of lattice Boltzmann models. Each column is a set of simulations under the designed cases. Each row is a set of simulations under selected lattice Boltzmann models including the theoretical solution.

图 4 各Taylor-Green涡设计计算工况下, 各格子Boltzmann模型计算误差随时间演变, 图中标号a, b, c, d分别对应表3各设计计算工况; 为方便表示, 图中横坐标单位为$ 0.173286795\; \rm{s} $
Figure4. The error evolution of each simulation. The panel a, b, c, d renders each designed case in Table 3 respectively. For the sake of rendering, the unit of the time axis is scaled as 0.173286795 s.

$ {\rm {Error}} = \frac{{\sqrt {\displaystyle\sum {{{\left( {u - {u_a}} \right)}^2} + {{\left( {v - {v_a}} \right)}^2}} } }}{{\sqrt {\displaystyle\sum {u_a^2 + v_a^2} } }}, $

式中$ u_a $和$ v_a $为理论结果.
从计算结果来看, 对于每个DnHm模型, 其数值精度随着$ u_{L{\rm B} ,0} $增大而变坏, 无论是计算结果等高线与理论解的相似程度(图3), 还是误差演化(图4)都体现了该趋势. 但与模型正值性一受到破坏计算即发散的D2H2模型不同的是, 高阶DnHm模型对于模型正值性破坏并未体现出间断式的数值表现变化. 对于单个计算工况, DnHm系列模型除不满足Galilean不变性的D2H2模型外, 计算精度随$ u_{\rm pos} $减小而变坏. 而D2H2模型精度则差于所有高阶模型, 包括正值区域更小的D2H4模型. 对比D2H2模型与其他高阶模型的数值表现, 模型的Galilean不变性对于模型的计算稳定性具有显著影响. 从这一角度而言, 矩精度对模型数值表现影响主要体现在模型是否满足Galilean不变性上. 对于满足Galilean不变性的模型, 矩精度影响不再显著, 模型数值表现主要受正值性影响.
分析模型正值区域内计算表现, 即$ \Delta u $大于0的数值模拟, 所有计算结果等高线与理论解均保持较好的符合, 而对比误差(图5), 计算结果未明显表现出模型正值区域与精度之间跨工况相关性. 具体而言, 在模型正值区域内, 具有较宽正值区域的模型并不对计算恒定保持优势精度, 以D2H3-2模型工况b为例, 其计算差于所有其他3阶及3阶以上DnHm模型工况a计算结果. 进一步而言, 在模型正值区域内, 其计算精度并不是恒定的, 计算工况越靠近正值区域上限, 模型精度越低. 因此模型正值区域优势主要体现在同一计算工况对比.
图 5 Taylor-Green涡设计计算工况下, $ \Delta u $值为正值的各格子Boltzmann模型数值模拟误差对比. 图中实线, 虚线和虚点线分别为工况a, b, c (见表3)结果, 离散模型则用颜色和符号标注. 为方便表示, 图中横坐标单位为0.173286795 s
Figure5. The error comparison of Taylor-Green vortex simulations with postive value of $ \Delta u $. The designed configurations of Taylor-Green vortex are denoted with line styles, in which the solid, dashed and dash-dotted line indicates case a, b and c in Table 3 respectively, meanwhile the lattice Boltzmann models are labeled with line colors and markers. For the sake of rendering, the unit of the time axis is scaled as 0.173286795 s.

在数值计算中, 独立的精度指标是非常实用的指示参量, 即对于任意计算工况, 只需保证该精度指标一致, 则计算精度一致. 鉴于正值区域的精度优势主要体现在同工况下模型之间的对比, 无法表征数值模拟的确切精度, 本文尝试探索其他参量作为精度指标的可能性. 在Taylor-Green涡算例中, 文章设计了$ \Delta u $值和$ \kappa $值来表征模型在对应计算工况下的正值性情况. 本文探索这两指标作为衡量计算精度的可能性. 文章首先检验$ \Delta u $值, 以D2H4模型工况a, b, c与D2H5模型工况b, c, d为例, 两两分别具有近似的$ \Delta u $值. 若格子Boltzmann方法的数值表现与$ \Delta u $值强相关, 则D2H4模型与D2H5模型对应算例精度将呈现明显的相关性. 从图6 (a)来看, 对比D2H3-2模型随工况变化而导致精度变化程度, 对应D2H4模型与D2H5模型结果之间的精度差异幅度依旧明显, 并不足以体现格子Boltzmann方法计算精度与$ \Delta u $值相关性. 另外取$ \Delta u $值在$ (-0.30\pm0.6) $附近计算结果 (图6 (b)), 同样对比D2H3-2模型随工况变化而导致精度变化程度, 选取的计算精度并未收敛于足够识别$ \Delta u $值影响的值域附近. 从上述分析可以看出, $ \Delta u $值并不具备衡量格子Boltzmann方法数值精度的能力. 对于$ \kappa $值, 本文选取了模型正值性已遭受破坏但初始$ \kappa $值尚未体现显著偏离的计算结果 (图7). 类似地, 对比D2H3-2模型随工况变化而导致精度变化程度, 误差对比显示, 选取算例并未显示显著的收敛性, 因此$ \kappa $值也不具备衡量格子Boltzmann方法数值表现的能力.
图 6 Taylor-Green涡数值模拟误差与$ \Delta u $值相关性分析　(a)对比了$ \Delta u $接近的D2H4模型a, b, c工况与D2H5模型b, c, d工况计算误差; (b)对比了$ \Delta u $在$ -0.30 $附近所有数值模拟计算误差. 所有图中均保留了D2H3-2模型在a, b, c, d工况下的计算误差作为参照. 曲线上标注值为该数值模拟的$ \Delta u $值. 图中a, b, c, d工况分别用实线, 虚线, 虚点线及点线标注, 而离散模型则用颜色和符号标注. 为方便表示, 图中横坐标单位为0.173286795 s
Figure6. The numerical performance of Taylor-Green vortex simulations vs the value of $ \Delta u $. Panel a plots the numerical errors of model D2H4 under case a, b, c and model D2H5 under case b, c, d, which possess close value of $ \Delta u $ respectively. Panel b plots the numerical errors of simulations with a value of $ \Delta u $ around $ -0.30 $. The numbers labeled on the curves are their values of $ \Delta u $. The simulation configurations are denoted with line style, in which solid, dashed, dash-dotted and dotted line indicates case a, b, c and d respectively, meanwhile the lattice Boltzmann models are labeled with line colors and markers. In all panels, the results of model D2H3-2 with four designed configurations are plotted for a reference. The designed configurations of Taylor-Green vortex are denoted with line styles, in which the solid, dashed, dash-dotted and dotted line indicates case a, b, c and d in Table 3 respectively, meanwhile the lattice Boltzmann models are labeled with line colors and markers. For the sake of rendering, the unit of the time axis is scaled as 0.173286795 s.

图 7 Taylor-Green涡数值模拟误差与初始$ \kappa $值相关性分析.　图中对比了在$ \Delta u $值为负情况下, 初始$ \kappa $值偏离$ 1 $幅值可忽略的算例计算误差, 包括模型D2H3-1工况c, 模型D2H3-2工况d, 模型D2H4工况b, 模型D2H5工况c, d和模型D2H6工况b, c. 图中均保留了D2H3-2模型在a, b, c, d工况下的计算误差作为参照. 曲线上标注值为该数值模拟的$ \Delta u $值. 图中a, b, c, d工况分别用实线, 虚线, 虚点线及点线标注, 而离散模型则用颜色和符号标注. 为方便表示, 图中横坐标单位为0.173286795 s
Figure7. The numerical performance of Taylor-Green vortex simulations vs initial $ \kappa $. The error evolutions of numerical simulations with negative values of $ \Delta u $ but negligible departure of initial $ \kappa $ from $ 1 $ are plotted, including model D2H3-1 under case c, model D2H3-2 under case d, model D2H4 under case b, model D2H5 under case c,d and model D2H6 under case b, c. The results of model D2H3-2 with four designed configurations are also plotted for a reference. The designed configurations of Taylor-Green vortex are denoted with line styles, in which the solid, dashed and dash-dotted line indicates case a, b and c in Table 3 respectively, meanwhile the lattice Boltzmann models are labeled with line colors and markers. For the sake of rendering, the unit of the time axis is scaled as 0.173286795 s.

对比经典的D2Q9模型, 在正值区域内, D2H2模型结果与之惊人的一致, 如工况a下两者的等高线图和误差演变曲线所示. 但在稳定性上, D2H2模型则略差于D2Q9模型. 而本文涉及的所有高阶模型(矩精度大于等于3)计算稳定性和精度都优于D2Q9模型. 从本文数值计算结果来看, D2H3-2模型具有最佳的计算稳定和精度, 同时对比其他高阶DnHm系列模型, 其计算量也具有明显优势. 因此从数值计算来说, D2H3-2模型值得进一步研究.

本文用Taylor-Green涡回避了高阶DnHm系列离散模型边界处理缺乏的现状, 数值验证了Ye等^[25]提出的离散模型的精度受到模型的正值性影响的推测. 文章选择了矩精度为{2—6}的最紧凑对称离散模型, 以及具有较宽正值区间的D2H3-2模型作为分析对象, 并与经典D2Q9模型进行比较. 根据涉及离散模型的正值区域, 文章通过调整时间步长设计了横跨各模型正值区域与非正值区域的四个计算工况. 分析各模型对四个计算工况数值模拟的结果, 文章得出了如下结论.
1) DnHm系列在正值区域内, 均能维持不错的计算结果, 但是并不能保证统一精度, 即随着计算工况靠近正值区域上限计算精度会随之下降. 在非正值区域计算时, 模型的表现受模型是否满足Galilean不变性影响. 具备Galilean不变性的模型, 其在非正值区域上也具有一定的计算稳定性, 并未出现间断式数值表现变化. 而不具备Galilean不变性的模型, 其计算精度在非正值区域存在间断式变化, 即非正值区域计算直接发散.
2)对于DnHm系列模型, 矩精度的影响主要体现在是否满足Galilean不变性. 对于不满足Galilean不变性模型, 除了其在非正值区域数值精度存在间断式变化外, 其正值区域内精度也差于正值区域更小的高阶模型. 而对于满足Galilean不变性的模型, 矩精度影响不再显著, 其数值表现主要取决于正值区域.
3)对于满足Galilean不变性的模型, 在相同计算工况下, 模型精度与模型正值区域正相关, 即模型正值区域越宽计算结果精度越好. 鉴于模型正值性仅适合分析相同计算工况模型的相对数值精度表现, 本文探索了$ \Delta u $值与$ \kappa $值用于衡量格子Boltzmann方法数值表现的可能性, 从分析结果来看, 两者并不具备表征能力.
4) D2H2模型与经典的D2Q9模型有着非常一致的精度表现, 尽管稳定性略有不如. 而高阶DnHm模型在稳定性和计算精度上都优于D2Q9模型. 特别是D2H3-2模型, 在所有本文涉及的DnHm模型中, 具有最佳的稳定性和精度, 是值得深入研究的模型.
平衡分布正值性对计算的影响可以从粒子分布需保持正值得到解释: 离散模型对于Maxwell分布全局正值性的破坏导致模型仅在有限范围有效, 该范围影响着模型的计算表现. 需要指出的是, 该有限范围并不等同于模型正值区域, 其原因是离散模型在正值区域内也仅为Maxwell分布的近似. 因此计算中, 模型正值区域内的精度并不一致. 而正值性与矩精度影响对比分析显示, 三阶Hermite多项式展开近似即能很好近似Maxwell分布函数形状. 另外, 本文对于离散模型的分析并不考虑模型间Ma数一致性问题, 即计算工况在不同模型下Ma数并不相同. 对于需要严格保持Ma数一致性的计算问题, 相关结论还需进一步验证分析, 本文在此不过多展开.
总结而言, 本文通过数值分析首次系统性梳理了离散模型正值性对计算影响, 并与矩精度进行对比分析. 文章证实了正值性对计算的影响, 确认通过提升模型正值性改善计算的可行性, 为pGHQ离散方案中可行方案筛选与未来模型进一步改进提供了方向.