顶盖驱动空腔流在石油开采、管道、运输等领域广泛存在, 深入研究其流动特性具有重要的理论意义和实际应用价值^[1]. 常见流动往往呈现湍流形态, 一般可以归结为高雷诺数下复杂流动^[2], 而这些问题都需要大量的计算, 其研究需要稳定有效的数值方法^[3].

LBM^[4,5] 将微观层次上的分子动力学模型与宏观层次上的流体动力学方程相结合, 将流体在时间、位置空间和速度空间上完全离散, 把流体视为一系列只有质量的粒子微团^[6], 这些流体粒子在给定的规则下相互碰撞和迁移, 整个过程遵循质量守恒和能量守恒. 与传统的数值计算方法相比, LBM 在模拟各种复杂流动问题方面具有程序简单、处理复杂边界灵活、图像清晰、计算局部性强和自然并行等优势^[7].

在 LBM 的许多模型中, 使用 Bhatnagar-Gross-Krook (BGK) 碰撞算子的单松弛时间 (SRT) 模型具有更简单的碰撞项, 是目前应用最广泛的一种 LB 模型^[8,9]. 然而, 该模型在模拟高雷诺数流时变得不稳定^[10]. 为此, D'Humieres 建立了MRT-LB 模型^[11], 在 MRT-LB 模型中, 碰撞过程在力矩空间中进行, 不同的弛豫时间对应不同的物理变量. MRT-LB 模型的数值稳定性和精度远远优于 SRT-LB 模型^[12], 而且该模型提供了更丰富的流场物理信息. 柴振华等^[13] 利用 Chapman-Enskog 多尺度展开方法、Mallwell 迭代、直接 Taylor 展开方法和递推方程方法从 LB 方程出发恢复出了 Navier-Stokes 方程. MRT-LB 模型已被广泛用于模拟流动, 尤其是高雷诺数流动^[14,15]. 景何仿等人^[16] 改进 MRT-LBM 模型, 并将其应用于刚性植被明渠流场的数值模拟. 张巧玲等人^[17] 应用 MRT-LBM 模拟了三维曲面腔的流动, 结果表明较少腔角的几何形状增强了稳定性. 由于传统的数值模拟方法不能稳定地模拟高雷诺数下的流动, 因此对于高雷诺数下的流动研究较少. 最近, 张凯等人^[18] 将 LES 思想引入到 MRT 中模拟了雷诺数在 500-10000 的方腔流. 已有研究表明随着雷诺数的增大腔内流动会进入到周期性状态^[19], 但对方腔内各流态的特征及跃迁临界值的研究较少. 然而, 当大单网格雷诺数很高时, 如果不增加额外的工作量, MRT-LB 模型就会变得不稳定^[20]. 因此, 在这项研究中, 将 MRT-LB 模型结合粘度抵消 (VC) 方法来解决这个问题.

Cheng 等人提出一种粘度抵消 (VC) 方法^[21,22], 其主要思想是在模拟具有较小粘度的实际流动中引入较大粘度所对应的人工粘度, 然后通过在方程中添加一个外力项来抵消该粘度, 从而稳定地模拟粘度较小的实际流动, 用于提高 SRT-LB 模型在高雷诺数流动模拟中的稳定性^[21]. 文献中已证实 MRT-LB 模型比 SRT-LB 模型具有更好的稳定性和准确性, 因此将 VC 方法引入 MRT-LB 模型可以获得更好的性能^[22]. 但 MRC-VC 方法可以模拟的最大雷诺数还需要进一步探究, 故本研究将利用 MRT-VC 方法研究雷诺数为 200-500000 的流动问题, 并将模拟结果与经典数值算例进行对比分析. 此外, 还重点研究了方腔流的几种流态以及不同流态之间的跃迁雷诺数临界值.

在不失去一般性的情况下, 我们考虑了具有 MRT 碰撞算子的二维九速 (D2Q9) 模型^[23]

其中, 离散速度 $\left \{ \boldsymbol{c}_{\alpha}:\alpha =0,1, \dots, 8 \right \}$ 定义为 $\boldsymbol{c}_{0} =\left ( 0,0 \right )$, $\boldsymbol{c}_{1} =-\boldsymbol{c}_{3} =\left ( 1,0 \right ) c$, $\boldsymbol{c}_{2} =-\boldsymbol{c}_{4} =\left ( 0,1 \right ) c$, $\boldsymbol{c}_{5} =-\boldsymbol{c}_{7} =\left ( 1,1 \right ) c$, $c=\Delta x/\Delta t$, $\Delta x$ 是格子空间步长, $\Delta t$ 是时间步长. $\hat{\boldsymbol{S}}$ 为非负对角矩阵. $\left\{\boldsymbol{m}_{\alpha}^{e q}(\mathbf{x}, t): \alpha =0,1,\cdots, 8\right \}$ 为平衡态分布函数. $\left \{ f_{\alpha }\left ( \boldsymbol{x},t \right ) : \alpha =0,1, \cdots,8 \right \}$ 是在 $\boldsymbol{x}$ 点 $t$ 时刻的离散分布函数. $\left \{\boldsymbol{m}_{\alpha}(\mathbf{x}, t):\alpha =0,1,\cdots, 8\right \}$ 为矩空间中的分布函数, $\left \{F_{\alpha}, \quad \alpha=0,1, \cdots, 8\right \}$ 为外力项 $\boldsymbol{F}$ 的强迫项. $\boldsymbol{M}$ 为变换矩阵

矩空间的平衡分布函数 $m^{e q}$ 为

松弛变换矩阵 $\hat{\boldsymbol{S}}$ 为

D2Q9 的松弛参数取 $s_{\rho}=s_{\varepsilon}=s_{\chi}=1. 0, s_{e}= s_{q}=1. 2$, $s_{v}=1/\tau$, $\tau$ 为松弛时间.

外力项 $F_{\alpha }$ 采用 Cheng 等人^[21] 提出的 LBE 力项, 可写成

其中 $g_{\alpha }\left ( \boldsymbol{x},t \right )$ 的计算公式为

上式中, $\boldsymbol{u}$ 是速度, $\omega_{\alpha}$ 是权重参数: $\omega_{0}=4/9,\omega_{1-4}=1/9,\omega_{5-8}=1/36$, $A$ 是流体连续性方程中的源项, 本文中 $A=0$. $F$ 是动量方程的外力项. 对于稳定流或变化缓慢的非稳定流, 为简单起见,下一时间步 $g_{\alpha}(\boldsymbol{x}+c_{\alpha }\Delta _{t},t+\Delta _{t})$ 由 $g_{\alpha}(\boldsymbol{x},t)$ 计算. 密度 $\rho$ 和速度 $\boldsymbol{u}$ 定义为

在高雷诺数的条件下, 流动的粘度较低, 数值模拟中的小耗散会导致计算不稳定. 纳维-斯托克斯 (N-S) 方程中的粘度项包含一个二阶导数, 它具有耗散作用,有助于维持数值稳定性. 在 N-S 方程的数值模拟中, 高阶截断误差在保持稳定性和精度方面发挥重要的作用. 广泛使用的 LB 模型在空间上是二阶的, 可以通过修改三阶或更高阶的截断误差项来提高稳定性, 而不会影响收敛速度. 另外, 在 LBM 中, 松弛时间不能小于 0.5, 即粘度必须为非负值. 然而, 在低粘度计算过程中, 由于数值误差的影响, 粘度误差可能为负, 这也是导致数值不稳定的重要原因之一^[21].

因此, 在 N-S 方程中增加一个额外的粘度项, 如式 (2.8) 所示, 可以抵消粘度误差, 并将其作为外部强制项进行抵消处理, 这种方法可以提高数值模拟的稳定性

在公式 (2.8) 的右侧, 第二项中的 $v+v_{c}$ 是模拟中指定的粘度 (应用粘度), $v_{c}$ 是需要抵消的粘度 (抵消粘度), 二者之差即为所需的模拟粘度 $v$ (目标粘度). 在粘度抵消法中, 外力 $\boldsymbol{F}_{i}$ 的计算公式为

其中 $S_{ij} =\frac{1}{2} \left ( \frac{\partial u_{i} }{\partial x_{j} } +\frac{\partial u_{j} }{\partial x_{i} } \right )$ 是应变速率张量. 因为 LB 模型在空间上是二阶的, 公式 (2.9) 的计算也应达到二阶或更高精度. $S_{ij}$ 可通过分布函数的非平衡部分计算得出

其中, $s_{v}$ 是与应用粘度相对应的松弛因子, 即

$S_{ij}$ 的偏导数可通过四阶精度的有限差分方案求得

需要注意的是, 对于边界附近的节点, $S_{ij}$ 的偏导数可通过外推法计算. 利用上述粘度抵消方法, 我们可以提高 MRT-LBM 模拟的稳定性. 在数值实验的基础上, 可以通过人为设置应用粘度来增强稳定性. 为了将结果误差限制在可接受的范围内, 应用的粘度不能太大.

MRT-VC 模型在模拟大雷诺数方腔流时有明显的优势, 为了更好的分析不同雷诺数下的方腔内流态, 本文使用 MRT-VC 模型模拟了不同雷诺数下的顶盖驱动方腔流.

本文计算区域如图1 所示, 其中矩形内部区域为流体区域. 顶部为驱动腔体内流体流动的平板, 图中的红点是流速监测点所在位置, 其坐标为 $\left ( 0. 2,0. 4 \right )$. 定义此流动的雷诺数为 Re$=U_p L_p/nu_p$, 其中$U_p=0. 5$m/s 为流体流动速度, $nu_p=1{\rm e}-6$ 为流体运动粘性系数, $L_p$ 为物理特征长度. 顶盖驱动速度 $U_l=0. 1$ 为平板运动速度, 下边界和两侧边界采用具有二阶精度的无滑移边界条件^{[24⇓⇓-27]}.

在数值模拟时将参数无量纲化为格子单位. 时间步长 $d_t=0. 001$, 节点总数 $N=100+1$, $nu_l$, $U_l$ 分别为无量纲化后的粘性系数和驱动速度. 无量纲化公式如下

本文设置 MRT-VC 模型中的粘度抵消系数为 0.5, 计算得到不同雷诺数下的抵消粘度和目标粘度如表1所示

为了检验 MRT-VC 能否准确稳定地模拟大雷诺数复杂涡流, 我们将其应用于顶盖驱动方腔流问题. 我们设置方腔上边界的速度 $U_{l} =0. 1$, 下边界和两侧边界采用无滑移条件. 首先, 我们使用 $256 \times 256$ 网格模拟 Re$=400,1000,2000,5000$ 的流动, 需要的目标粘度 $v$ 分别为 $0. 0640,0. 0256,0. 0128,0. 0051$, 模拟结果以流线展示为图2, 可以发现随着雷诺数的增大, 主涡涡心逐渐向腔体的集合中心移动, 这与文献 [18] 的结果一致. 当雷诺数为 5430 时, SRT 模型仅可以模拟腔内中心附近主涡 (初级涡) 的流线型结构, 不能模拟二次涡, 当雷诺数进增加到 7000 时, SRT 模型无法再模拟初级涡和次级涡^[18]. 然而, 从图2 可以看出, 应用粘度 ($v+v_{c}$) 为 0.0055 的 MRT-VC 模型可以精确模拟中心附近的初级涡及腔壁附近的次级涡流, 当雷诺数为 7000 时, 腔内右下角附近生成三级旋涡, 这与 Ghia 等人^[25] 的结果一致, 因此可以得出结论, MRT-VC 模型稳定更好, 可以用于模拟高雷诺数下的湍流. 表2 对 Re=400、1000、2000 和 5000 流线下的主涡和副涡中心进行了定量比较, 其中 (a), (b), (c), (d) 分别代表本文、文献 [1,24,25] 中的结果, 也说明了 MRT-VC 模型的准确性.

为了进一步验证 MRT-VC 模型是否增强了数值稳定性, 我们将流动雷诺数提高到 Re=20000, 目标粘度 $v$ 设置为 0.0013, 其他参数不变. 在这些条件下, 流动变得湍急, 出现了更多的附属涡旋, 这与已有结果^[26]一致. 涡旋在各自的周期内演化和运动, 没有出现不稳定现象, 表明 MRT-VC 模型具有更强的数值稳定性. 此外利用 MRT-VC 模型, 我们可以模拟 Re= 60000、100000 和 200000 的流动, 如图3 所示. 随着雷诺数的增大, 初级涡形状不规则, 附属涡旋数目逐渐增加.

矩形区域的网格数设置为 $100 \times 100$. 当雷诺数较小时, 流场一般处于稳态流流动状态, 但在 MRT-VC 模型中具体稳态流场的雷诺数的范围是多少, 需要继续探究. 图4 为方腔内稳定流动状态时流场内流函数等值线的分布情况. 从图中可以看出, 方腔内主涡位于右上部分区域. 随着雷诺数的增加, 主涡逐渐向腔体的几何中心靠近, 右下角、左下角、左上角的次涡尺寸不断变大.

为了探究流速随时间的变化, 在计算过程中, 我们在流场中选取坐标为 $(0. 2,0. 4)$ 的点来监测该点速度随时间的变演化过程, 图5 为该点处 $x$ 方向和 $y$ 方向速度分量随时间演化过程及相空间内速度图谱, 可以看出, 速度在某一时刻之后逐渐稳定, 空间内速度轨迹逐向一个点靠近, 称之为速度稳定点. 当 Re $=200$ 时, 速度分量以震荡的方式逐渐趋于稳定, 速度轨迹类似于半径不断减小的圆, 最后达到稳定点的速度为 $(-0. 011848,-0. 00291887)$, 当 Re $=10000$ 时, 速度分量逐渐减小并趋于平衡, 速度轨迹逐渐向左下角线性减小, 最后达到稳定点的速度为 $(-0. 038219,-0. 013007)$.

图6 为不同雷诺数下的轴向中心流速分布. 由图可见, 稳态流的轴向流速呈线性分布 (除边界附近外), 雷诺数越大, 最大流速与最小流速的差值越大.

随着流动雷诺数的继续增大, 方腔内流动将不再呈现出稳定状态, 进入到周期性的变化过程中. 从图7 可以看出周期性流态的方腔内次涡尺寸较雷诺数最大的稳态流没有明显的增大, 但在右下角、左下角、左上角逐步生成了三级次涡.

图8 说明进入周期态的方腔流轴向速度分布的线性特征衰减. 图9 可以看出当 Re=12500 时, 监测点横向和纵向速度均呈周期性均匀变化, 速度轨迹为一条简单闭曲线, 当 Re=45000 时, 监测点的速度分布存在一种特殊的周期性, 随着时间的推进每个周期内速度分布的形状相同, 但速度峰值逐渐降低, 相空间速度轨迹呈螺旋状向左下角行进.

当雷诺数增大到 Re=50000 时, 横向速度仍然存在一定的周期性, 如图10 所示, 横向速度仍然随时间周期性降低, 但每个周期内的速度演化呈锯齿状, 纵向速度没有明显的变化规律, 先横向速度进入混沌流动状态, 将此流动状态称为不完全混沌状态流场. 图11 展示了不完全混沌状态下的方腔流流函数等值线分布. 可以看出, 在该状态下, 主涡发生形变, 雷诺数越大, 形变越严重. 轴向流速严重偏离线性分布, 如图12 所示.

当雷诺数为 Re=50000 时, 主涡产生轻微形变, 横向速度随时间震荡式下降, 速度轨迹不规则地螺旋前进; 当雷诺数为 Re=95000 时, 主涡产生严重的形变, 有向右上角分离的趋势, 横向速度和相空间图谱的演化明显呈现出多个周期.

随着流动雷诺数增大到 Re=100000, 主涡开始分离, 二级次涡发生形变, 方腔内流动将进入到混沌流动状态, 混沌状态方腔流的流函数等值线分布如图13. 可以看出, 混沌态方腔流内小涡的数量随着雷诺数的增加而增多, $X$ 方向和 $Y$ 方向的速度分量呈震荡分布, 雷诺数越大, 震荡越明显, 如图14.

使用 MRT-VC 模型能模拟的最大雷诺数 Re=500000, 此时腔内的流动情况如图13, 可以看出, 此时流场内流动非常复杂, 主涡分离出若干小涡, 右下角二级次涡消失, 左上角流动有明显的不稳定. 相空间图谱如图15 所示, 不难发现, 中线流速震荡剧烈, 相空间内速度轨迹变化呈现明显的混沌状态.

综上, 在同一参数下, 方腔内流态从稳态流到周期态流的雷诺数临界值在 10000-12500 之间, 从周期态流到半周期-半稳态流的雷诺数临界值在 45000-50000 之间, 从半周期-半混沌态流到混沌流的雷诺数临界值在 95000-100000 之间.

本文通过在 MRT-LBM 模型中引入粘度抵消方法, 建立了 D2Q9 MRT-VC 模型. 我们使用该模型对二维顶盖驱动方腔流进行了数值模拟研究, 计算了雷诺数在 200-500000 范围内流场的流动情况. 从仿真结果中, 可以得到一些重要结论

$\bullet$ 本文将改进的 MRT-VC 模型的仿真结果与其他参考文献中的模拟结果进行了对比分析, 验证了本文所建立模型的准确性, 并模拟了大雷诺数条件下的流动情况. 模拟结果表明 MRT-VC 模型可显著提高计算的稳定性, 适用于模拟研究大雷诺数流动问题.

$\bullet$ 我们探究了 MRT-VC 模型可以模拟的最大雷诺数, 发现该模型可模拟的最大雷诺数达到 500000. 在模拟较大雷诺数的流场时, 模拟结果表明: 主涡涡心逐渐向方腔的几何中心移动, 然后发生形变, 最后分离出若干个小涡, 二级次涡尺寸逐渐变大, 三级次涡逐渐形成, 最后左下角的次涡尺寸缩小且右下角的次涡发生形变并消失.

$\bullet$ 当雷诺数不超过 10000 时, 方腔内的流场为稳定态流场. 稳定态流场的流线光滑, 主涡呈圆形, 流速分量随时间增大逐渐趋向于一个固定的值, 速度轨迹逐渐缩小为一个点, 轴向流速呈线性分布.

$\bullet$ 当雷诺数在 12500-45000 时, 流场为周期态流场. 周期态流场内的主涡中心位于腔体的几何中心, 三级次涡生成, 速度分量随时间周期性变化, 速度轨迹为一条简单闭曲线或螺旋线.

$\bullet$ 当雷诺数在 50000-95000 时, 流场为不完全混沌态流场. 不完全混沌态流场内主涡发生形变, 有分离的趋势, 速度分量在每个周期内分布不均匀或者呈现出多个周期, 轴向速度脱离线性分布.

$\bullet$ 当雷诺数大于 100000 时, 流场为混沌态流场. 混沌态流场主涡分离出诺干个小涡, 横向速度分量和速度轨迹呈明显的混沌状态, 轴向流速震荡剧烈.

$\bullet$ 我们发现从稳定流到周期流的临界跃迁雷诺数在 10000-12500 之间, 周期流到不完全混沌态流的临界跃迁雷诺数在 45000-50000 之间, 不完全混沌态流到混沌流的临界跃迁雷诺数在 95000-100000 之间.

$\bullet$ 当雷诺数超过 100000 时, 使用 MRT-VC 模型仿真方腔流出现局部物理震荡现象, 方腔左上角的流线出现一个有限的粗糙区域, 后续将改变粘度抵消比例系数、细化左上角网格、调整边界条件来进一步改善此模型.