时间步长对格子玻尔兹曼法模拟室内气流精度的影响

摘要

基于格子玻尔兹曼法的大涡模拟（LBM-LES）是湍流模拟的新方法，但不恰当的时间步长${\delta }_t$可能会影响其计算精度。首先理论总结了${\delta }_t$可能对LBM-LES湍流模拟造成的影响，阐明过大的${\delta }_t$会导致速度场产生压缩性误差，而过小的${\delta }_t$会导致超松弛碰撞产生速度场的数值振荡。其次，通过对等温室内气流案例进行LBM-LES模拟，定量讨论了${\delta }_t$引起的压缩性误差和数值振荡问题。结果表明，${\delta }_t$较大时流场密度变化剧烈，且格子玻尔兹曼单位的马赫数（M）超过0.3的区域中速度场产生了明显的压缩性误差。同时，过小的${\delta }_t$导致平均及脉动风速均产生了数值振荡，这在网格分辨率较高时尤为明显。建议模拟时在确保${\delta }_t$足够小以满足最大风速区域的M<0.3的基础上，尽量增大${\delta }_t$以防止产生数值振荡。

关键词

风环境; 格子玻尔兹曼法; 时间步长; 压缩性误差; 超松弛碰撞; 大涡模拟

近年，基于格子玻尔兹曼法的大涡模拟（lattice Boltzmann method-based large-eddy simulation， LBM-LES）开始应用于建筑 ^［1-3］和城市风环境^［4］模拟。与当前风环境主流的有限体积法（finite volume method， FVM）在宏观尺度上求解物理量不同，LBM用虚拟的、包含有限种速度模式的微观分布函数表示流体粒子的集合，并通过分布函数的碰撞和迁移来模拟流体运动^［5］。与基于FVM的大涡模拟 （FVM-LES，即风环境模拟的主流LES方法）^［6-7］相比，LBM-LES算法简单，边界条件易于实现，且在LES计算中无需求解压力泊松方程^［8］，计算速度更快，在复杂湍流风环境模拟中具有较大潜力^［9-10］。

LBM的控制方程是格子玻尔兹曼方程，其中的关键项是碰撞算子。碰撞算子的形式决定了待求解流体的性质。BGK（Bhatnagar-Gross-Krook）近似模型^［11］是最常用的碰撞算子，能以简单的形式获得足够精度的模拟结果。既往研究表明，BGK算子在模拟环境风问题，尤其是室内气流中取得了较好的成果。Elhadidi等^［1］利用含BGK算子的LBM模拟了较粗网格条件下室内空间的气流分布并与FVM进行了比较。Han等^［2］则系统总结了使用LBM-LES模拟室内空气流动时不同计算条件对模拟精度的影响，并与FVM-LES进行了详细比较。上述研究表明含BGK算子的LBM-LES可有效模拟室内空气流动，并可取得与FVM-LES类似的模拟结果。

同时理论分析表明，利用含BGK算子的LBM-LES可推导出低马赫数（Mach number， M）流体的连续性方程和纳维‒斯托克斯（Navier-Stokes， N-S）方程，但推导出的N-S方程与风环境模拟中常用的方程形式有所差异^［12-14］，从而产生模拟误差。该误差与模拟时间步长${\delta }_t$的设定紧密相关，不恰当的${\delta }_t$设置可能导致计算结果出现压缩性误差或数值振荡。目前，在应用LBM-LES模拟风环境问题时，${\delta }_t$的设置引起的误差问题尚未引起足够重视和充分讨论。为了在应用LBM-LES模拟风环境时对如何设置${\delta }_t$提供依据，本文梳理了既往研究，从理论方面系统总结了${\delta }_t$的设置可能引起的压缩性误差或数值振荡，并以室内空气流动案例为例，定量讨论了${\delta }_t$造成压缩性误差或数值振荡的程度。

1 含BGK的LBM中与时间步长相关的误差理论分析与总结

1.1　含BGK的格子玻尔兹曼方程及LBM-LES方法的理论回顾

本节首先对含BGK碰撞算子的格子玻尔兹曼方程及LBM-LES方法的理论进行简要回顾，以便后文对误差进行理论分析。该方程如式（1）所示。

$$\begin{array}{l}{f}_a\left(\mathit{r}+{\delta }_t{\mathit{e}}_a,t+{\delta }_t\right)-f_a\left(\mathit{r},t\right)=\\ -\frac{\mathrm{1}}{\tau }\left[f_a\left(\mathit{r},t\right)-f_a^{\mathrm{e}\mathrm{q}}\left(\mathit{r},t\right)\right]\end{array}$$

（1）

式中：$f_a$为a方向分布函数；${\mathit{e}}_a$为a方向上$f_a$的离散速度；$\mathit{r}$和$t$分别为$f_a$所在位置向量和时间；$f_a^{\mathrm{e}\mathrm{q}}$为$f_a$的平衡函数；${\delta }_t$为离散时间步长；$\tau$为$f_a$的松弛时间。

格子玻尔兹曼法方程在微观层面描述了流体粒子的分布函数随时间发展的演化过程。当分布函数确定后，流体速度$\mathit{u}$、密度$\rho$及压力$p$等宏观物理量可通过式（2）求得。其中，$e_{\mathrm{s}}$为格子声速，在三维问题中的值为$\mathrm{1}/\sqrt[]{\mathrm{3}}$，其他参数含义同前。

$$\rho =\sum _a{f}_a\left(\mathit{r},t\right)$$

（2a）

$$\mathit{u}=\frac{\mathrm{1}}{\rho }\sum _a{\mathit{e}}_a{f}_a\left(\mathit{r},t\right)$$

（2b）

$$p=\rho e_{\mathrm{s}}^{\mathrm{2}}$$

（2c）

在面对高雷诺数Re的湍流问题时，可基于LBM开展LES计算（LBM-LES）。根据LES理论，流体的总粘性${\nu }_{\mathrm{t}\mathrm{o}\mathrm{t}}$由分子粘性$\nu$及亚格子粘性${\nu }_{\mathrm{s}\mathrm{g}\mathrm{s}}$共同构成（即${\nu }_{\mathrm{t}\mathrm{o}\mathrm{t}}=\nu +{\nu }_{\mathrm{s}\mathrm{g}\mathrm{s}}$）^［15］。同时，基于LBM理论，流体的总粘性${\nu }_{\mathrm{t}\mathrm{o}\mathrm{t}}$与总松弛时间${\tau }_{\mathrm{t}\mathrm{o}\mathrm{t}}$存在如式（3）所示的关系：

$${\nu }_{\mathrm{t}\mathrm{o}\mathrm{t}}={e_{\mathrm{s}}}^{\mathrm{2}}({\tau }_{\mathrm{t}\mathrm{o}\mathrm{t}}-\mathrm{0.5})$$

（3）

基于式（3）可用总松弛时间${\tau }_{\mathrm{t}\mathrm{o}\mathrm{t}}$替换式（1）中的松弛时间$\tau$以开展LBM-LES计算。与传统FVM-LES相同，只需采取合适的LES亚网格模型计算亚网格粘性${\nu }_{\mathrm{s}\mathrm{g}\mathrm{s}}$，即可开展LBM-LES计算。

1.2　实际物理量到格子玻尔兹曼单位的转换

在式（1）中，以真实物理量（通常包含量纲）所度量的流体原型问题首先被映射到格子玻尔兹曼单位的物理量（通常是量纲一的）进行模拟，模拟完成后再将其映射回真实物理量以输出结果。故模拟第一步是确定适当的转换参数。在无外力等温流体问题中，LBM主要关注流体粘性$\nu$、速度$\mathit{u}$、压力$p$及位置$\mathit{r}$参数，在进行上述物理量的转换时只需网格分辨率${\delta }_x$和时间步长${\delta }_t$两个转换参数。各物理量转换关系如式（4）所示，其中上标$\mathrm{p}\mathrm{h}$和$\mathrm{l}\mathrm{b}$分别表示真实物理量及其对应的格子波尔兹曼单位物理量。

$${\mathit{r}}^{\mathrm{l}\mathrm{b}}={\mathit{r}}^{\mathrm{p}\mathrm{h}}\frac{\mathrm{1}}{{\delta }_x}$$

（4a）

$${\mathit{u}}^{\mathrm{l}\mathrm{b}}={\mathit{u}}^{\mathrm{p}\mathrm{h}}\frac{{\delta }_t}{{\delta }_x}$$

（4b）

$$p^{\mathrm{l}\mathrm{b}}=p^{\mathrm{p}\mathrm{h}}{\left(\frac{{\delta }_t}{{\delta }_x}\right)}^{\mathrm{2}}$$

（4c）

$${\nu }^{\mathrm{l}\mathrm{b}}={\nu }^{\mathrm{p}\mathrm{h}}\frac{{\delta }_t}{{{\delta }_x}^{\mathrm{2}}}$$

（4d）

其中，网格分辨率${\delta }_x$通常根据湍流复杂度、模拟所需精度及计算量共同确定，与传统FVM-LES基本相同。需要注意的是，LBM中采用的是均匀正立方体网格，无法如FVM-LES一样在局部复杂湍流区域加密网格，故应注意使用测试网格独立性等方式来确定网格尺寸。这在既往研究中已得到多次验证^［2，4］。而时间步长${\delta }_t$的确定则与FVM-LES有较大差异。过大或过小的${\delta }_t$都可能导致明显的精度误差，这将是本文接下来的讨论重点。

1.3　压缩性误差

利用BGK算子，可从格子玻尔兹曼方程中推导出形如式（5）的N-S方程^［14-15］：

$$\frac{\partial \left(\rho \mathit{u}\right)}{\partial t}+\nabla \cdot \left(\rho \mathit{u}\mathit{u}\right)=-\nabla p+\nabla \cdot e_{\mathrm{s}}^{\mathrm{2}}\left(\tau -\frac{\mathrm{1}}{\mathrm{2}}\right)\rho \left[\nabla \mathit{u}+{\left(\nabla \mathit{u}\right)}^{\mathrm{T}}\right]+\mathrm{O}\left({ϵ}^{\mathrm{2}}\right)+\mathrm{O}\left(M^{\mathrm{3}}\right)$$

（5）

式（5）中：$\mathrm{O}\left({ϵ}^{\mathrm{2}}\right)$和$\mathrm{O}\left(M^{\mathrm{3}}\right)$分别为与${ϵ}^{\mathrm{2}}$和$M^{\mathrm{3}}$相关的高阶省略项；$ϵ$为努森数（Knudsen number），与马赫数M成正相关，即$\mathrm{O}\left({ϵ}^{\mathrm{2}}\right)~\mathrm{O}\left(M^{\mathrm{2}}\right)$ ^［12］。式（5）相比标准的N-S方程多了与M相关的附加项。其中M是以格子波尔兹曼单位定义的，如式（6）所示：

$$M=\frac{\left|{\mathit{u}}^{\mathrm{l}\mathrm{b}}\right|}{e_{\mathrm{s}}}=\frac{\left|{\mathit{u}}^{\mathrm{p}\mathrm{h}}\right|{\delta }_t}{e_{\mathrm{s}}{\delta }_x}$$

（6）

式中：$\left|{\mathit{u}}^{\mathrm{p}\mathrm{h}}\right|$为局部物理速度${\mathit{u}}^{\mathrm{p}\mathrm{h}}$的大小。

式（5）同时表明，推导出的N-S方程以可压缩形式存在（即无法消去各项中的密度$\rho$），意味着LBM-LES在处理不可压缩流体问题时本质上是一种伪可压缩方法，并会产生所谓的“压缩性误差”^［12］。虽然严格来说不可压的流体并不存在，但在处理诸如建筑与城市风环境等低流速问题时，FVM-LES通常采用不可压缩的N-S方程。由此可见LBM-LES与FVM-LES在计算时针对是否可压缩的处理方法上有较大区别。既往研究^［12，15］表明，LBM的压缩性误差包括与密度梯度相关的误差，以及上述与M相关附加项导致的误差。

式（6）表明，即使在局部风速与网格分辨率${\delta }_x$一定时，较大的时间步长${\delta }_t$会增大M，从而增大使式（5）中与M及$ϵ$相关项的值（即压缩性误差）。故在模拟不可压缩湍流问题时，与传统FVM-LES相比，不恰当的${\delta }_t$可能导致M增大，从而使得模拟结果产生压缩性误差。Skordos^［16］曾尝试用LBM模拟层流状态下的二维泰勒涡旋流和剪切流，发现涡旋流的模拟值和解析值之间的误差随着M的减小而减小，并最终实现稳定收敛。而在剪切流中，随着M的减小，模拟误差呈先减小而后增大的趋势。Reider等人^［12］从理论上推导了压缩性误差，并证实在Re=100且周期为2π的衰减泰勒涡流模拟结果准确性随着M的降低而提高。

应当注意，BGK算子中的压缩性误差与FVM-LES中库朗数C的不正确设置引起的误差并非同一概念，尽管C也是由${\delta }_x$和${\delta }_t$间的取值关系造成。在FVM-LES中，在处理低M不可压缩流体时，通常建议选择适当的时间步长${\delta }_t$以将C控制在小于1（即$C=\left|{\mathit{u}}^{\mathrm{p}\mathrm{h}}\right|{\delta }_t^{\mathrm{p}\mathrm{h}}/{\delta }_x^{\mathrm{p}\mathrm{h}}<\mathrm{1}$），否则将造成结果误差甚至模拟发散。而在LBM-LES中，保证模拟稳定性的一个必要条件是M<1，即$\left|{\mathit{u}}^{\mathrm{l}\mathrm{b}}\right|<\mathrm{1}/\sqrt[]{\mathrm{3}}\approx \mathrm{0.577}$，否则模拟将直接发散。故若要保证LBM-LES模拟正常稳定进行，则必有$C=\left|{\mathit{u}}^{\mathrm{l}\mathrm{b}}\right|{\delta }_t^{\mathrm{l}\mathrm{b}}/{\delta }_x^{\mathrm{l}\mathrm{b}}<\mathrm{1}$始终成立（因为LBM规定了${\delta }_t^{\mathrm{l}\mathrm{b}}=\mathrm{1}$和${\delta }_x^{\mathrm{l}\mathrm{b}}=\mathrm{1}$）。

1.4　超松弛与数值振荡

从式（1）中可明显看出BGK算子的本质是表现分布函数$f_a\left(\mathit{r},t\right)$以一定的速率向平衡分布状态$f_a^{\mathrm{e}\mathrm{q}}\left(\mathit{r},t\right)$的演化过程，即松弛过程。该式可改写为如式（7）的形式：

$$\begin{array}{l}{f}_a\left(\mathit{r}+{\delta }_t{\mathit{e}}_a,t+{\delta }_t\right)=\\ \left(\mathrm{1}-\frac{\mathrm{1}}{\tau }\right)f_a\left(\mathit{r},t\right)+\frac{\mathrm{1}}{\tau }{f}_a^{\mathrm{e}\mathrm{q}}\left(\mathit{r},t\right)\end{array}$$

（7）

根据$\left(\mathrm{1}-\frac{\mathrm{1}}{\tau }\right)$的取值为正、负或零，$f_a$或缓慢接近$f_a^{\mathrm{e}\mathrm{q}}$，或立刻达到$f_a^{\mathrm{e}\mathrm{q}}$，或直接超过$f_a^{\mathrm{e}\mathrm{q}}$。BGK算子可导致分布函数有如下三种演化形态^［17］：

（1） 当$\mathrm{1}-\frac{\mathrm{1}}{\tau }>\mathrm{0}$，即$\tau >\mathrm{1}$时，$f_a$以固定速率向$f_a^{\mathrm{e}\mathrm{q}}$逐渐演化，称为亚松弛（under relaxation）；

（2） 当$\mathrm{1}-\frac{\mathrm{1}}{\tau }=\mathrm{0}$，即$\tau =\mathrm{1}$时，$f_a$只需一个时间步长即达到$f_a^{\mathrm{e}\mathrm{q}}$，称为全松弛（full relaxation）；

（3） 当$\mathrm{1}-\frac{\mathrm{1}}{\tau }<\mathrm{0}$，即$\frac{\mathrm{1}}{\mathrm{2}}\le \tau <\mathrm{1}$时，$f_a$直接超过$f_a^{\mathrm{e}\mathrm{q}}$，称为超松弛（over relaxation）。

应当注意，$\tau$不可小于$\frac{\mathrm{1}}{\mathrm{2}}$，因为根据式（3），流体的粘性不可为负。Krüger 等人^［17］研究了初始条件为$f_{\mathrm{0}}/f_{\mathrm{0}}^{\mathrm{e}\mathrm{q}}=\mathrm{1.1}$、且$f_{\mathrm{0}}^{\mathrm{e}\mathrm{q}}$为恒定值条件下的BGK算子，得到了如图1所示的$f_{\mathrm{0}}$与$f_{\mathrm{0}}^{\mathrm{e}\mathrm{q}}$的关系，对应了上述的三种形态模式。

图 1  BGK算子中的亚松弛、全松弛及超松弛算例（重绘自Krüger等 ^［17］）

Fig. 1  Example of under, full, and over relaxation in BGK collision (reproduced from Krüger et al ^[17])

图1表明理想的碰撞过程是亚松弛或全松弛，即$f_a^{}$平滑地或直接向$f_a^{\mathrm{e}\mathrm{q}}$演化。在实际模拟中，全松弛难以达到，因为不可能经过一步就完成模拟。而在超松弛中$f_a$将围绕$f_a^{\mathrm{e}\mathrm{q}}$振动并以指数幅度衰减，最后达到$f_a^{\mathrm{e}\mathrm{q}}$。但$\tau$过小可能导致振动过于剧烈，$f_a$无法达到$f_a^{\mathrm{e}\mathrm{q}}$，最终得到错误结果。

综合式（3）、式（4）可得到如式（8）的关系：

$${\tau }_{\mathrm{t}\mathrm{o}\mathrm{t}}=\frac{{\nu }_{\mathrm{t}\mathrm{o}\mathrm{t}}^{\mathrm{l}\mathrm{b}}}{{e_{\mathrm{s}}}^{\mathrm{2}}}+\mathrm{0.5}=\frac{{\nu }_{\mathrm{t}\mathrm{o}\mathrm{t}}^{\mathrm{p}\mathrm{h}}{\delta }_t}{{e_{\mathrm{s}}}^{\mathrm{2}}{{\delta }_x}^{\mathrm{2}}}+\mathrm{0.5}$$

（8）

式（8）表明，由于流体粘性${\nu }_{\mathrm{t}\mathrm{o}\mathrm{t}}^{\mathrm{p}\mathrm{h}}$通常为定值，故当确定网格分辨率${\delta }_x$后，${\tau }_{\mathrm{t}\mathrm{o}\mathrm{t}}$的大小与${\delta }_t$正相关。${\delta }_t$的设定影响${\tau }_{\mathrm{t}\mathrm{o}\mathrm{t}}$的大小，从而决定了碰撞过程的松弛模式。理想状态下$\tau$应不小于1，据式（8）可知需要${\nu }_{\mathrm{t}\mathrm{o}\mathrm{t}}^{\mathrm{l}\mathrm{b}}\ge \frac{\mathrm{1}}{\mathrm{6}}$或${\nu }_{\mathrm{t}\mathrm{o}\mathrm{t}}^{\mathrm{p}\mathrm{h}}\ge \frac{{{\delta }_x}^{\mathrm{2}}}{\mathrm{6}{\delta }_t}$。然而这在风环境模拟中较难满足。由于空气动粘性系数极小（数量级为${\mathrm{10}}^{-\mathrm{5}}{\mathrm{m}}^{\mathrm{2}}·{\mathrm{s}}^{-\mathrm{1}}$），即使在LES计算中加上亚网格粘性${\nu }_{\mathrm{s}\mathrm{g}\mathrm{s}}$也很难大于$\frac{{{\delta }_x}^{\mathrm{2}}}{\mathrm{6}{\delta }_t}$。例如，在风环境模拟中，由于库朗数和计算量的双重限制，很难在设置${\delta }_x=$1 m的同时使得${\delta }_t={\mathrm{10}}^{-\mathrm{5}}$ s；也无法使用${\delta }_x={\mathrm{10}}^{-\mathrm{3}}$ m来匹配${\delta }_t={\mathrm{10}}^{-\mathrm{1}}$s。尽管在风速较小的局部区域可能满足${\nu }_{\mathrm{t}\mathrm{o}\mathrm{t}}^{\mathrm{p}\mathrm{h}}\ge \frac{{{\delta }_x}^{\mathrm{2}}}{\mathrm{6}{\delta }_t}$，但建筑或城市尺度的风环境中，大部分区域（特别是所关注的区域）主要以高雷诺数湍流为主，故风环境模拟中BGK算子主要表现为超松弛模式。通常这种超松弛模式下的振动类似湍流脉动，然而不恰当的网格分辨率${\delta }_x$和时间步长${\delta }_t$之间的关系会导致类似图1中$\tau =\mathrm{0.51}$所示的剧烈振动，最终导致模拟结果产生数值振荡。尤其是在相同的网格分辨率${\delta }_x$下，过小的${\delta }_t$会导致${\nu }^{\mathrm{l}\mathrm{b}}$过小，从而导致严重的数值振荡。这与传统的FVM-LES具有极大的不同。

综合上述分析可知，在LBM-LES进行风环境模拟时，当${\delta }_x$确定后，过大的${\delta }_t$可能导致与M相关的项产生较大的压缩性误差，而${\delta }_t$过小则使松弛时间$\tau$减小，可能造成松弛碰撞算子产生数值振荡。这是LBM-LES相比传统FVM-LES在模拟设置上的一个重要区别。FVM-LES中，只要离散时间步长满足C<1则不会对模拟结果产生显著的影响。既往研究已经讨论了层流状态下二维流动中LBM-LES的压缩性误差^［12］，而以风环境模拟为代表的湍流状态下不同时间步长对模拟结果的影响尚未充分讨论。下节将以室内气流为例定量讨论这一问题。该案例边界条件相对简单、纯粹、易于控制，便于进行定量研究。

2 等温室内气流模拟案例

本文采用国际能源机构推荐的标准等温室内气流案例（IEA-Annex 20^［18］），对含BGK算子的LBM-LES进行压缩性误差研究。该案例形状及取样位置如图2所示，房间特征参数为$L/H=\mathrm{3}$， $h/H=\mathrm{0.056}$ ，$r/H=\mathrm{0.16}$。其中L， W和H分别代表房间长度、进深和高度，且$H=\mathrm{3.0}$ m。h和r分别代表气流入口及出口高度。由房间高度及气流入口速度定义的Re≈89 000。模拟参数详表1。本文中，含标准Smagorinsky亚网格模型^［19］及BGK碰撞算子的LBM-LES应用于本模拟。Han等^［2］的既往研究已经表明该亚网格模型和BGK算子能较好地适应室内气流模拟，取得满意的模拟精度。

图 2  等温室内气流案例的几何形状及取样点分布（修改自IEA-Annex 20 ^［18］）

Fig. 2  Geometry and sampling points distribution of isothermal indoor airflow case (modified from IEA-Annex 20 ^[18])

表 1  模拟参数及相关边界条件

Tab. 1  Simulation parameters and boundary conditions

参数项目	参数值
亚网格模型	标准Smagorinsky模型 （Cs = 0.12 ［20］）
房间尺寸	9.0 m× 3.0 m × 3.0 m （L × W × H）
LES模拟流场物理时间	预备模拟区间： 18 min，采样区间： 6 min（换气率 0.172 min-1）
入口边界条件	$U_{\mathrm{i}\mathrm{n}}=0.455$ m∙s-1，无流入湍流
出流边界条件	压力梯度0
固体壁面边界条件	Bounce-back边界（无滑移）

之前Han等^［2］ 已经讨论了LBM-LES的网格分辨率对模拟精度的影响，本文基于其研究结果仅选取${\delta }_x=\mathrm{1}/\mathrm{75}H$及$\mathrm{1}/\mathrm{150}H$两种网格，讨论不同时间步长${\delta }_t$对模拟精度的影响。具体工况设置如表2所示。工况名称规定为“XaTb”，表示网格分辨率及时间步长分别为${\delta }_x=H/a$及${\delta }_t=\mathrm{1}/b$ s。

表 2  工况设置

Tab. 2  Case settings

工况名称	${\delta }_x$/ m	${\delta }_t$ / s	${\nu }^{\mathrm{l}\mathrm{b}}$		工况名称	${\delta }_x$/ m	${\delta }_t$ / s	${\nu }^{\mathrm{l}\mathrm{b}}$
X75T50	0.04 m （1/75 H）	1/50	$1.82\times {10}^{-4}$		X150T200	0.02 m （1/150 H）	1/200	$9.10\times {10}^{-5}$
X75T100		1/100	$9.10\times {10}^{-5}$		X150T400		1/400	$4.55\times {10}^{-5}$
X75T200		1/200	$4.55\times {10}^{-5}$		X150T800		1/800	$2.28\times {10}^{-5}$
X75T400		1/400	$2.28\times {10}^{-5}$		X150T1600		1/1 600	$1.14\times {10}^{-5}$
X75T800		1/800	$1.14\times {10}^{-5}$		X150T3200		1/3 200	$5.70\times {10}^{-6}$
X75T1600		1/1 600	$5.70\times {10}^{-6}$

3 模拟结果与讨论

3.1　平均与湍流脉动风速结果及精度分析

图3显示了$\overline{u}$（时间平均风速的x方向分量）和$\sqrt[]{\overline{u{\mathrm{\text{'}}}^{\mathrm{2}}}}$（基于时间平均的脉动风速标准差的x方向分量）的模拟结果。所有结果均用入口风速$U_{\mathrm{i}\mathrm{n}}$进行量纲归一化。图中添加了Nielsen等^［21］的实验数据用于验证模拟的准确性。

图 3  部分区域$\overline{\mathit{u}}$与$\sqrt[]{\overline{\mathit{u}{\mathbf{\text{'}}}^{\mathrm{2}}}}$的量纲一化模拟结果与实验结果对比

Fig. 3  Comparison of simulation and experiment results of normalized $\overline{u}$ and $\sqrt[]{\overline{u{\mathrm{\text{'}}}^{\mathrm{2}}}}$ in some regions

除精度最低的X75T50，几乎所有工况模拟结果均能再现$\overline{u}$和$\sqrt[]{\overline{u{\mathrm{\text{'}}}^{\mathrm{2}}}}$的分布趋势，且模拟精度有随着${\delta }_t$的减小而提高的趋势。X75T800和X75T1600中，$\overline{u}$和$\sqrt[]{\overline{u{\mathrm{\text{'}}}^{\mathrm{2}}}}$均出现了轻微的空间数值振荡。而在X150工况组中，X150T200的$\overline{u}$和$\sqrt[]{\overline{u{\mathrm{\text{'}}}^{\mathrm{2}}}}$都达到了最佳精度。随着${\delta }_t$的减小，精度并未提高，反而出现了明显的数值振荡，从而降低模拟精度。

本文采用式（9）所定义的L2误差范数 ^［17］${\epsilon }_q$定量评估模拟精度。式中$q_{\mathrm{E}\mathrm{X}\mathrm{P}}\left(\mathit{r}\right)$和$q_{\mathrm{L}\mathrm{B}\mathrm{M}}\left(\mathit{r}\right)$分别代表实验和模拟中位置$\mathit{r}$处的物理量q。L2误差范数考虑了Nielson实验数据的所有点。${\epsilon }_q$越小代表模拟与实验之间的误差越小，从而模拟精度越高。

$${\epsilon }_q=\sqrt[]{\frac{{\sum }_{\mathit{r}}{\left(q_{\mathrm{L}\mathrm{B}\mathrm{M}}\left(\mathit{r}\right)-q_{\mathrm{E}\mathrm{X}\mathrm{P}}\left(\mathit{r}\right)\right)}^{\mathrm{2}}}{{\sum }_{\mathit{r}}{q_{\mathrm{E}\mathrm{X}\mathrm{P}}}^{\mathrm{2}}\left(\mathit{r}\right)}}$$

（9）

图4显示了不同${\delta }_t$时所有工况的L2误差范数变化曲线。随着${\delta }_t$从1/50 s减小到1/200 s，X75工况组中$\overline{u}$和$\sqrt[]{\overline{u{\mathrm{\text{'}}}^{\mathrm{2}}}}$的误差均减小，随后则显著增大。可以推测，${\delta }_t$从1/50 s减小到1/200 s时的精度提升可能是由于压缩性误差的减小；而${\delta }_t<\mathrm{1}/\mathrm{200}$ s时精度降低应该是由于超松弛导致，因为在这些工况中观察到了$\overline{u}$和$\sqrt[]{\overline{u{\mathrm{\text{'}}}^{\mathrm{2}}}}$的数值振荡。对于X150工况组，${\delta }_t=\mathrm{1}/\mathrm{200}$ s时，模拟精度最高，而后随着${\delta }_t$的减小$\sqrt[]{\overline{u{\mathrm{\text{'}}}^{\mathrm{2}}}}$的模拟精度迅速衰减，这也应归因于超松弛引起的振荡。3.2和3.3节将深入讨论这些推测。

图 4  不同${\mathit{\delta }}_{\mathit{t}}$的模拟精度L2误差范数变化曲线

Fig. 4  L2 error norms of different ${\delta }_t$ values

3.2　压缩性误差的讨论

X75T50、X75T100和X75T200这三个工况的模拟精度有较大差异，但其风速模拟结果并未显示出明显的数值振荡，且三者的网格设置完全一致仅${\delta }_t$不同。这表明它们的精度差异极有可能是由于不同${\delta }_t$导致的压缩性误差所引起，于是本节选择这三个工况分析压缩性误差。

图5显示了三个工况中各区域的时间平均密度$\overline{\rho }$与初始值${\rho }_{\mathrm{0}}$相比的相对偏差$\left(\overline{\rho }-{\rho }_{\mathrm{0}}\right)/{\rho }_{\mathrm{0}}$。如1.3节所述，LBM-LES是一种伪可压缩模拟方法，即使在模拟同一个不可压缩问题时，不同的${\delta }_t$设置亦会导致密度变化显著。X75T50的${\delta }_t$较大，导致$\overline{\rho }$明显偏离了初始值，尤以入口附近区域更为明显。

图 5  部分工况的中央垂直平面量纲一化时间平均密度偏差$\left(\overline{\rho }-{\rho }_{\mathrm{0}}\right)/{\rho }_{\mathrm{0}}$分布

Fig.5  Deviation of time-averaged density $\left(\overline{\rho }-{\rho }_{\mathrm{0}}\right)/{\rho }_{\mathrm{0}}$ at central vertical section in some cases

图6显示了所有工况中整个模拟域的空间平均密度相对于初始值的偏差。该偏差反映了LBM-LES计算中密度的压缩程度。随着${\delta }_t$减小一半，空间平均密度的差异几乎呈指数衰减。当密度偏差小于0.5 %后逐渐达到稳定。此时的密度接近初始值，可忽略压缩性的影响。

图 6  所有工况中全空间平均密度与初始值的相对偏差变化

Fig. 6  Deviation between spatial-averaged density and initial value of all cases

垂直截面上量纲一化$\overline{u}$及M的计算结果如图7所示。在X75T50中，来自入口的气流有远离天花板向下的趋势，导致该区域的模拟结果精度较差。据图7b所示，该区域M大于该工况其他区域及其他工况的M，并超过了0.3。这与Krüger等^［17］的研究一致。该研究建议LBM-LES模拟场中的由${\mathit{u}}^{\mathrm{l}\mathrm{b}}$定义的M不应大于0.3，否则将产生较显著的压缩性误差。同时，图5显示该区域密度变化剧烈。这再次表明较大的密度梯度将会导致显著的压缩性误差。随着M的降低，入口处的气流趋于水平，表明压缩性误差可通过降低M得到一定的补偿。X75T50的其他区域或其他工况中M均小于0.3，故可忽略$\overline{u}$的压缩性误差。同时，X150工况组中M<0.3始终成立，表明X150工况组的压缩性误差均不明显，因而X150工况组中模拟精度并未随着${\delta }_t$的降低而提高。

图 7  部分工况中垂直截面上量纲一化$\overline{\mathit{u}}$及M分布

Fig. 7  Vertical distribution of normalized $\overline{u}$ and M of some cases

以上讨论证实，在使用LBM-LES求解室内湍流时，过大的M会导致速度场产生明显的压缩性误差。通过调整${\delta }_t$可减小M，以补偿误差。为了减少压缩性误差造成的影响，应尽量保证流场中最大风速区域的M < 0.3，即${\delta }_t\le \frac{\sqrt[]{\mathrm{3}}{\delta }_x}{\mathrm{10}\left|{\mathit{u}}^{\mathrm{p}\mathrm{h}}\right|}$。值得注意的是，M是由${\mathit{u}}^{\mathrm{l}\mathrm{b}}$定义的格子玻尔兹曼单位的参数。即使模拟问题原型相同，也可以通过使用不同的${\delta }_t$改变M，这与物理场中由真实速度定义的M不同。

3.3　超松弛导致的数值振荡讨论

3.1节图3显示，在X75T400、X75T800、X75T1600工况和X150工况组中的大多数工况中均观察到明显数值振荡，这可能由于超松弛碰撞模式造成，本节将对此进行分析。图8显示了所有工况的$f_{\mathrm{0}}$在点$\left(x,y,z\right)=(H,H/\mathrm{2},H/\mathrm{2})$处当流场达到充分发展后某个1 s周期（第1 080~1 081 s）内的松弛状况。选择点$\left(x,y,z\right)=(H,H/\mathrm{2},H/\mathrm{2})$是因为在该点处X75T50、X75T100和X75T200中没有明显的振荡，而在其他工况中发生振荡，即存在一个振荡发生的临界状况。图中$f_{\mathrm{0}}(t)/f_{\mathrm{0}}^{\mathrm{e}\mathrm{q}}(t)$表明某一$t$时刻$f_{\mathrm{0}}$与$f_{\mathrm{0}}^{\mathrm{e}\mathrm{q}}$间的大小关系

图 8  所有工况在$\left(\mathit{x},\mathit{y},\mathit{z}\right)=(\mathit{H},-\mathit{H}/\mathrm{2},\mathit{H}/\mathrm{2})$位置处1 s内${\mathit{f}}_{\mathrm{0}}$的超松弛现象

Fig. 8  Over relaxation phenomenon in one second at $\left(x,y,z\right)=(H,-H/\mathrm{2},H/\mathrm{2})$ of all cases

在所有工况中，$f_{\mathrm{0}}$围绕$f_{\mathrm{0}}^{\mathrm{e}\mathrm{q}}$来回摆动，表明所有工况都是超松弛碰撞模式，而不是理想的亚松弛模式，此时$f_{\mathrm{0}}$并不向$f_{\mathrm{0}}^{\mathrm{e}\mathrm{q}}$呈指数衰减。这一结果证实了在湍流中，超松弛碰撞模式比亚松弛更常见。在X75工况组中，$f_{\mathrm{0}}$摆动的“频率”随${\delta }_t$的减小而增大。同时摆动“幅度”随${\delta }_t$减小而减小。从X75T50到X75T200，摆动过程似乎没有形成数值振荡，而应该是由湍流脉动导致。然而在X75T400、X75T800和X75T1600中，摆动过程演化为可见的高频振动，与速度场的模拟结果发生数值振荡的状况一致。类似地，在X150组中，当${\delta }_t$降低到一定程度时，分布函数形成了高频振动，最终形成了速度场发生数值振荡，这在X150T1600至X150T3200之间尤为明显。

由此可见，当${\delta }_t$减小到一定程度时，超松弛碰撞模式所对应的湍流脉动会最终演化成分布函数的高频振动，并最终导致宏观速度场的数值振荡。然而，超松弛碰撞模式何时演化为高频振动则较为复杂，其与局部湍流的流动模式、网格尺寸、流体性质及碰撞算子等都有很大关系，并非只与${\delta }_t$线性相关，故较难判断数值振荡的临界${\delta }_t$。一般建议在消除压缩性误差的前提下尽量增大${\delta }_t$以避免数值振荡。

4 结论

本文分析并总结了采用含BGK碰撞算子的LBM-LES模拟风环境问题时，时间步长${\delta }_t$对模拟精度的影响，并以室内气流模拟案例对其进行了定量讨论。主要结论如下：

（1）LBM-LES是一种伪可压缩方法，在处理不可压缩问题时模拟域中的密度在模拟过程中会发生变化。过大的${\delta }_t$会使得密度变化较大，导致速度场产生压缩性误差。较小的${\delta }_t$可以减小压缩性误差。

（2）在模拟湍流时，BGK碰撞算子通常表现为超松弛碰撞模式，即分布函数表现为一定程度的摆动。过小的${\delta }_t$会导致摆动会演化成高频振动，最终使得速度场发生数值振荡。该现象在网格分辨率相对较高时更容易产生。

（3）在确定网格分辨率（如网格独立性测试）后，${\delta }_t$应足够小以满足最大风速区域M < 0.3（即${\delta }_t\le \frac{\sqrt[]{\mathrm{3}}{\delta }_x}{\mathrm{10}\left|{\mathit{u}}^{\mathrm{p}\mathrm{h}}\right|}$）以减小压缩性误差。在此基础上尽量采用较大的${\delta }_t$以防止数值振荡的发生。

本文仅粗略确定了${\delta }_t$的取值上限，今后的工作将着重于如何确定${\delta }_t$的合理范围，并建立${\delta }_t$与其他物理量之间的定量关系。此外，对应于风环境中高Re问题的模拟，通常采用比BGK更为复杂、鲁棒性更高的碰撞算子（如MRT、cumulant LBM等），${\delta }_t$的变化对这些碰撞算子的影响也应予以进一步研究。

作者贡献声明

韩梦涛：制定研究目标及内容，数值模拟，数据分析，论文撰写及修订。